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ABSTRACT 

We present numerical simulations of the spectral evolution and emission of radio components in 
relativistic jets. We compute jet models by means of a relativistic hydrodynamics code. We have 
developed an algorithm (SPEV) for the transport of a population of non-thermal electrons including 
radiative losses. For large values of the ratio of gas pressure to magnetic field energy density, ~ 
6 X 10^, quiescent jet models show substantial spectral evolution, with observational consequences only 
above radio frequencies. Larger values of the magnetic field {a^ ~ 6 x 10^), such that synchrotron 
losses are moderately important at radio frequencies, present a larger ratio of shocked-to-unshocked 
regions brightness than the models without radiative losses, despite the fact that they correspond to 
the same underlying hydrodynamic structure. We also show that jets with a positive photon spectral 
index result if the lower limit 7min of the non-thermal particle energy distribution is large enough. 
A temporary increase of the Lorcntz factor at the jet inlet produces a traveling perturbation that 
appears in the synthetic maps as a super luminal component. We show that trailing components can 
be originated not only in pressure matched jets, but also in over-pressured ones, where the existence 
of recollimation shocks does not allow for a direct identification of such features as Kelvin-Helmholtz 
modes, and its observational imprint depends on the observing frequency. If the magnetic field is large 
(ttg ^ 6 X 10^), the spectral index in the rarefaction trailing the traveling perturbation does not change 
much with respect to the same model without any hydrodynamic perturbation. If the synchrotron 
losses are considered the spectral index displays a smaller value than in the corresponding region of 
the quiescent jet model. 

Subject headings: galaxies: jets - hydrodynamics - radiation mechanisms: nonthermal - relativity 
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1. INTRODUCTION 

Relativistic jets are routinely observed emerging from 
active galactic nuclei and microquasars. and presumably 
they arc behind the phenomenology detected in gamma- 
ray bursts. It is a broadly recognized fact that the ob- 
served VLBI radio-maps of parsec-scale jets are not a 
direct map of the physical state (density, pressure, veloc- 
ity, magnetic field) of the emitting plasma. The emission 
structure is greatly modified by the fact that a distant 
(Earth) observer detects the radiation emitted from a jet 
which moves at relativistic speed and forms a certain an- 
gle with respect to the line of sight. Time delays between 
different emitting regions, Doppler boosting and light 
aberration shape decisively the observed aspect of every 
time-dependent process in the jet. The observed patterns 
are also influenced by the travel path of the emitted ra- 
diation towards the observer since Faraday rotation and, 
most importantly, opacity modulate total intensity and 
polarization radio maps. Finally, there arc other effects 
that can be very important for shaping VLBI observa- 
tions which do not unambiguously depend on the hydro- 
dynamic jet structure, namely, radiative losses, particle 
acceleration at shocks, pair formation, etc. In this work 
we try to account for some of these elements by means 
of numerical simulations. 

The basis for currently accepted interpretation of 
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the phenomenology of r elati vistic jets wa s set by 
iBlandford fc Konigll (IOTI and|Konigl| (fTMll) . A num- 
ber of analytic works have settled the basic under- 
standing that accounts for the non-thermal synchrotron 
and inverse Compto n emission of extragalactic jets 
(e.g., lMarschei1ll98C( ). as well as the spectral evolution 
of superluminal compone n ts in parsec-scale jets (e.g. , 
Blandford fc McKeel Il976l: iHiighes. Aller & Alleil fl985l : 



Marscher. Gear fc Travid 119921 : iMarscher fc Geailll985l) . 
Assuming kinematic jet models, the numerical implemen- 
tation of these analytic results enables one to extensively 
test the most critical theoretical assumptions by compar- 
ison wi th the observed p h enom e nology both for steady 
(e.g., iDalv fc Marscheil [T988t iHughes. AUct fc AUed 
1989a[ 119911: iGomez. Alberdi fc j^c ^aida [l993l [l994l : 



Gomez, et al.lll994D and unsteadv lets fe.g.. lJonesl ll988^. 
Basically, the aforementioned numerical implementation 
consists on integrating the synchrotron transfer equa- 
tions assuming that radiation originated from an ideal- 
ized jet model and accounting for all the effects men- 
tioned in the previous paragraph. 

The advent of multidimensional relativistic (magneto- 
) hydrodynamic numerical codes has allowed to replace 
the previously used kinematic, steady jet models by mul- 
tidimensional, time-dependent hydrodyn amic models o f 
thermal plasmas ( for a rev iew see, e.g., 'G6 mezl[200^ . 
The works of Gomez, et al.l fl995„ .1997.), (hereafter G95 
and G 97, respectivelv) iDuncan. Hughes fc OppermanI 
()1996() or iKomissarov fc Falld (|1996[ ) compute the svn- 
chrotron emission of relativistic hydrodynamic jet models 
with suitable algorithms that account for a number of rel- 
ativistic effects (e.g., Doppler boosting, light aberration, 
time delays, etc.). Their models assume that there ex- 
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ists a proportionality between the number and the energy 
density of the hydrodynamie (thermal) plasma and the 
corresponding number and energy density of the emit- 
ting population of non-thermal or supra-thermal parti- 
cles. These authors assumed that the magnetic field was 
dynamically negligible, that the emitted radiation had 
no back-reaction on the dynamics, and that that syn- 
chrotron losses were negligible. All these assumptions are 
very reasonable for VLBI jets at radio observing frequen- 
cies if the jet magnetic field is sufficiently weak. Consis- 
tent with their assumptions, the former papers included 
neither a consistent spectral evolution of the non-thermal 
particle (NTP) population, nor the proper particle and 
energy transport along the jet. 

The spectral evolution of NTPs and its transport in 
cl assical jets and radiogalax i es have been carried out 
bvl Jones. Rvu fc Engell (ll999l ).[Mic ono. et al.l (|1999[) and 
iTregillis. Jones fc Rvul ()200lD . In these works a coupled 
evolution of a non-relativistic plasma along with a pop- 
ulation of NTPs has been used to asses either the sig- 
natures of diffusive shock acceleration in radio galaxies 
(| Jones. Rvu fc Engel|[T999l : ITregillis. Jones fc Rvull2001h 
or the observational imprint of the non-linear saturation 
of Kelvin-Hel mholtz (KH) modes developed by a per- 
turbed beam (jMicono. et al.|[T999[ ). ICasse fc MarcowithI 
()2003f l have also developed a scheme to perform mul- 
tidimensional Newtonian magneto-hydrodynamical sim- 
ulations coupled with stochastic differential equations 
adapted to test particle acceleration and transport in 
kilo-parsec scale jets. Dealing with the spectral evolution 
of NTPs is relevant in view of the multiband observations 
of extragalactic jets where, a significant aging of the emit- 
ting particles s eems to be present at optical to X-ray fre- 
quen cies fM8 7.lHeinz fc Begelmanlll997HMarshall. et all 
120021: Cen A, iKraft. Forman. Jones fc Murravll200i[r 

This paper builds upon the lines opened by G95 and 
G97. G95 concentrated on the emission properties from 
steady relativistic jets, focusing on the role played by 
the external medium in determining the jet opening an- 
gle and presence of standing shocks. G97 used a similar 
numerical procedure to study the ejection, structure, and 
evolution of supcrluminal components throu gh variations 



in the ejection velocity at the jet inlet. lAgudo et al.l 
(|200lh discussed in detail how a single hydrodynamie 
perturbation triggers pinch body modes in a relativis- 
tic, axisymmctric beam which result in observable su- 
pcrluminal feat ures trailing t he ma in supcrluminal com- 
poncnt. Fina l lv. |Alov. et al.l ()2003l ) extended the work of 
lAgudo et al.l (|2001f ) to three-dimensional, helically per- 
turbed beams. Here, we combine multidimensional rela- 
tivistic models of compact jets with a new algorithm to 
compute the spectral evolution of supra-thermal parti- 
cles evolving in its bosom, i.e., including their radiative 
losses, and their relevance for the emission and the spec- 
tral study of relativistic jets. 

This work is composed of two parts. In the first part, 
wc present a new numerical scheme to evolve popula- 
tions of relativistic electrons in relativistic hydrodynam- 
ical flows including radiative losses (§ [3]) . For the pur- 
pose of calibration the new method, our work is based 
upon the same axisymmctric, relativistic, hydrodynamie 
jet models as employed in G97. Using the same jet pa- 
rameters allows us to quantify the relevance of including 
radiative losses and, along the way, to compare the emis- 



Table 1 

Set of models used in this work. Values in the table 

REFER to the JET NOZZLE. THE FIRST COLUMN LISTS THE 

model names. the second and third columns give 
the jet-to-external-medium pressure ratio, and the 

comoving magnetic field at the jet nozzle, 
respectively. the last column lists the values of 
the ratio of gas pressure to magnetic pressure. 
Additional models not including radiative losses, 

but computed with the spev method, will be 
denoted with a suffix "-nl" . likewise, models with 
the same parameters as the ones listed here, but 
computed with the am method, will be denoted 
with a suffix "-am". 
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6 X 10"' 


PM-H 
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6 X 10^ 


OP-L 


1.5 


0.03 


6 X 10"' 


OP-H 


1.5 


0.30 


6 X 10^ 



sion properties of parsec-scalc jets computed according 
to two different methods: (1) the new method presented 
in this paper and (2) the method presented in G95 and 
G97; to which wc will refer, for simplicity, as Adiabatic 
Method (AM). In the second part of the paper, we apply 
the new method to quantify the relevance of radiative 
losses in the evolution of both quiescent and dynamical 
jet models. We will show (§[5]) the regimes in which both 
approaches yield similar synthetic total intensity radio 
maps and when synchrotron losses modify substantially 
the results. We also show which are the key parameters 
to trigger a substantial NTP aging and, therefore, to sig- 
nificantly change the appearance of the radio maps cor- 
responding to the same underlying, quiescent jet models. 
The spectral evolution of a hydrodynamie perturbation 
travelling downstream the jet, will be discussed in Sect. [71 
Finally, we discuss our main results and conclusions in 
Sect. El 

2. HYDRODYNAMIC MODELS 

Two quiescent, relativistic, axisymmctric jet models 
constitute our basic hydrodynamie set up (see Tab. [1]). 
They correspond to the same pressure-matched (PM), 
and over-pressured (OP) models of G97. The models 
were comput ed in cylindrical sym metry with the code 
RGENESIS (|Mimica. et al.ll2004D . The computational 
domain spans (10i?b x 200i?b) in the (r x z)-plane (i?h 
is the beam cross-sectional radius at the injection posi- 
tion). A uniform resolution of 8 numerical cclls/i?h is 
used. The code module that integrates the relativistic 
hydrodynamics equations is a conservative, Eulerian im- 
plementation of a Godunov-type scheme with high-order 
spati al and temporal accuracy (based on the GEN ESIS 
code: lAlov. et al.lll999l: lAlov. Pons fc Ibaiiezlll999[) . We 
follow the same nomenclature as G97 where quantities 
affected by subscripts a, b and p refer to variables of 
the atmosphere, of the beam at the injection nozzle and 
of the perturbation (§ 12. ip . respectively. The jet ma- 
terial is represented by a diffuse {ph/pa = 10"'^; p be- 
ing the rest-mass density), relativistic (Lorentz factor 
Ft ~ 4) ideal gas of adiabatic exponent 4/3, with a 
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Mach number Mi, — 1.69. At the injection position, 
model PM has a pressure Pb = Pa, while model OP 
has Pb = l.SPa-Pressure in the atmosphere decays with 
distance z according to P{z) = Pa/[l + (z/zc)^'^]^'^^, 
where Zc = 6QRb- With such an atmospheric profile both 
jet models display a paraboloid shape, which introduces 
a small, distance-dependent, jet opening angle which is 
compatible with observations of parse-scale jets. At a 
distance of 200i?b, the opening angles for the models PM 
and OP are 0.29° and 0.43°, respectively. 

Pressure equilibrium in the atmosphere is ensured by 
including adequate counter-balancing, numerical source 
terms. However, despite the fact that the initial model 
is very close to equilibrium, small numerical imbalance 
of forces triggers a transient evolution that decays into a 
final quasi-steady state after roughly 2 — 5 longitudinal 
grid light-crossing times. We treat these quiescent states 
as initial models. Model PM yields an adiabatically- 
expanding, smooth beam. Model OP develops a collec- 
tion of cross shocks in the beam, whose spacing increases 
with the distance from the jet basis. 

2.1. Injection and Evolution of Hydrodynamic 
Perturbations 

Variations in the injection velocity (Lorcntz factor) 
have been suggest ed as a way to generate internal shocks 
in relativistic jets (|Reeslll978l) . We set up a traveling per- 
turbation in the jet as a sudden increase of the Lorentz 
factor at the jet nozzle (from Fb = 4 to Pp = 10) for a 
short period of time (0.75i?b/c; c being the light speed). 
Since the injected perturbation is the same as in G97, 
its evolution is identical to the one these authors showed 
and, thus, we provide a brief overview here. The per- 
turbation develops two Ricmann fa ns emerging from its 
leading and rear edges (see, e.g., iMimica. et al.l [20051 : 
IMimica. Alov fc Miillei1l2007l l. In front of the perturba- 
tion a shock-contact discontinuity-shock structure {SCS) 
forms, while the rear edge is trailed by a rarefaction- 
contact discontinuity-rarefaction (TZCTZ) fan. In the lead- 
ing shocked region the beam expands radially owed to the 
pressure increase with respect to the atmosphere. In the 
trailing rarefied volume the beam shrinks radially on ac- 
count of the smaller pressure in the beam than in the 
external medium. This excites the generation of pinch 
body modes in the beam that seem to trail the main hy- 
drody namic perturbation as pointed out bv lAgudo et al.l 
()2001h . Also the component itself splits in, at least, two 
parts when the forward moving rarefaction leaving the 
rear edge of the component merges with the reverse shock 
traveling backwards (in the component rest frame) that 
leaves from the forward edge of t he hydrodynamic per- 
turbation (as in I Alov. et al.ll2003l ). 

3. SPEV: A NEW ALGORITHM TO FOLLOW 
NON-THERMAL PARTICLE EVOLUTION 

The spectral evolution (SPEV) routines are a set of 
methods developed to follow the evolution of NTPs in 
the phase space. Here we assume that the radiative 
losses at radio frequencies are negligible with respect 
to the total thermal energy of the jet at every point 
in the jet. Thus, radiation back reaction onto the hy- 
drodynamic evolution is neglected. Certainly, such an 
ansatz is invalid at shorter wavelengths (optical. X-rays), 
where radiative losses shape the observed spectra (see. 



e.g., IMimica. et all (pOOl for X -ray-synchrotron blazar 
models that include the radiation back-reaction onto the 
component dynamics). 

The 7-dimensional space formed by the particle mo- 
menta, particle positions and time is split into two parts. 
For the spatial part of the phase space, we assume that 
NTPs do not diffuse in the hydrodynamic (thermal) 
plasma. Thereby, the spatial evolution of the NTPs is 
governed by the velocity field of the underlying fluid, and 
it implies that the NTP comoving frame is the same as 
the thermal fluid comoving frame. Assuming a negligible 
diffusion of NTPs is a sound approximation in most parts 
of our hydrodynamic models since the electron diffusion 
lengths are much smaller than the dynamical lengths 
in smooth flows (sec, e.g.. iTregillis. Jones fc Rvull2001l : 
[Miniati 200ii ). Obviously, the assumption is not fulfilled 
wherever diffusive acceleration of NTPs takes place (e.g., 
at shocks or at the jet lateral boundaries). Nevertheless, 
there exists a strong mismatch between the scales rele- 
vant to dynamical and diffusive transport processes for 
NTPs of relevance to synchrotron radio-to- X-ray emis- 
sions within relativistic jets. The mismatch ensures that 
even in macroscopic, non-smooth regions such as the 
cross shocks in the beam of model OP, the assumption 
we have made suffices to provide a good qualitative de- 
scription of the NTP population dynamics. 

Consistent with the hydrodynamic discretization, we 
assume that the velocity field is uniform inside each nu- 
merical cell (equal to the average of the velocity inside 
such cell). In practice, a number of Lagrangian parti- 
cles are introduced through the jet nozzle, each evolv- 
ing the same NTP distribution but being spatially trans- 
ported according to the local fluid conditions. We em- 
phasize that these Lagrangian particles are used here 
for the solely purpose of representing the spatial evo- 
lution (i.e., the trajectories) of ensembles of NTPs. We 
integrate the trajectories of such particles using a con- 
ventional time-explicit, adaptive-step-size, fourth order 
Runge-Kutta (RK) integrator. 

3.1. Particle Evolution in the momentum space 

In order to derive the equations governing the time evo- 
lution of charged NTPs i n the momentum space we follow 
closely the appr oach of iMiralles. van Riper fc Lattinie^ 
^99M (see also IWebbI I1985L or lKirklll994f) . We start 
by considering the Boltzmann equation that obeys the 
ensemble averaged distribution function / of the NTPs, 
each with a rest-mass toq ; 

where / is a function of the coordinates and the com- 
ponents of the particle 4-momentum with respect to 
the coordinate basis e(Q) . The F^^ are the usual Christof- 
fel symbols and the right hand side represents the colli- 
sion term, r being the particle proper time. 

Equation [1] can be written in terms of the particle 
4-momcntum components with respect to the comov- 
ing or matter frame instead of the components with re- 
spect to the coordinate basis. The comoving tetrad e^a) 
{a = 0, 1,2,3), is formed by four vectors, one of which 
(e(o)) is the four velocity of the matter and the following 
orthonormality relation is fulflUcd 

e(a) ■ = riab. 
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where rjab is the Minkowski metric {rjoo — — !)■ Wc ex- 
phcitly point out that the components of tensor quanti- 
ties with respect to the coordinate and tetrad basis are 
annotated with Greek and Latin indices, respectively. 
The transformation between the basis e(„) and e(^a) is 
given by the matrix e" and its inverse matrix e^, 

e(a) =^ e^e(„), e(„) = e'^e(^a) (2) 

In terms of the comoving basis, the Bohzmann equa- 
tion is 



P 



dl_ 



^ bcP 



dp'' 



(3) 



coll 



The connection coefficients in the tetrad frame F^^ obey 
the foUowing relations 



pa pi 



(4) 



where the comma and the semicolon stand for partial 
and covariant derivatives, respectively. 

We introduce the two first moments of the distribution 
function by the equations 

^2 



(5) 



(6) 



where p^ = — mgC^ is the square of the NTP three- 
momentum measured by the comoving observer. The 
solid-angle integrations are performed over all parti- 
cle momentum directions. The number of NTPs per unit 
volume with modulus of their three-momentum between 
p and p+dp for an observer comoving with the matter 
is n'^{p)dp. Further integration of the above moments 
n" and t"'' over p. dp, gives the hydrodynamical mo- 
ments. 

In order to obtain the continuity equation for NTPs, 
we multiply the Boltzmann equation ([3]) by (p^/p*^), and 
integrate oy er Q to yield (for the details see App. A of 
IWebbl[l985h . 



d 



(9p V P 



-i-ab 



coll 

The next step is to formulate the continuity equation 
in the diffusion approximation. Such approximation im- 
plies that the scattering of NTPs by hydromagnetic tur- 
bulence results in a quasi isotropic distribution function 
in the scattering (comoving) frame. Thus, it is assumed 
that the distribution function of the NTPs can be ex- 
pressed as the sum of two terms, / ~ f'--^^ + f^^^fl, where 
jio) ^ ^"(1) and fl is the unit vector in the direction of 
the momentum of the particle. With such an assump- 
tion, we obtain that 



p2 S'^ 



which also leads to 



(8) 



(9) 

Plugging the approximations ^ and ^ into Eq. ([7]) 
and neglecting the terms coming from the anisotropy of 



the distribution function, i.e., the terms arising from f^-^"^ 
we obtain 



0;q' 




Equation (jTU]) is valid for any general metric g^,y. How- 
ever, in the present work we are only interested in ob- 
taining the transport equation for NTPs in the special 
relativistic regime. To restrict Eq. (fTU]) to such a regime 
we take a fiat metric, (jfj,i/ — ^/ii/- Thereby, the tetrad 
and the coordinate frame basis are related by a simple 
Lorcntz transformation, i.e., 

p°-F 



= 5., + (F-l) 



(«,j = 1,2,3), 



where v"^ {i = 1,2,3) are the spatial components of the 
velocity of matter, which is equal to the hydrodynami- 
cal velocity of the NTPs, since we make the assumption 
that NTPs do not diffuse in the hydrodynamic plasma. 
The hydrodynamical Lorcntz factor of the plasma is de- 
noted by F = 1 / ■\/(l — v'^Vi). With this transformation 
we obtain 



dt 



' 'dt ' 



dTv' 
dx' 



= e 



(11) 



(12) 



O being the expansion of the underlying thermal fluid, 
which is related to p by 

D\np 



Dt 



(13) 



Plugging Eqs. ((II])-(IT3l) into Eq. ^ and using the 
definition of the Lagrangian derivative with respect to 
the proper time of the comoving observer 

D „ / a , d 



yields 



Dt 



d 



dt 



-pe 



L»T " dp\'i 
which can be cast in the form 

D\nn° p^ainn° 2^ 1 

— R + = — 

Dt 3 op 3 iv' 



dx^ 



dn- 



, (14) 



coll 



dn 



P° \^^/ coll 

The collision term contains the interaction between 
NTPs and matter, radiative losses due to synchrotron 
processes, etc. Let us consider first the interaction with 
matter. In this case, the collisions can be assumed to 
be isotropic in the comoving frame and elastic. In such 
a case, and consistently to the previous approximation, 
the collision term in Eq. (jisp vanishes and we can find 
a solution for the homogeneous differential equation by 
considering 

D In n'^ p 9 In n° 
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as the derivative of In nP along the following curve in the 
plane (t,p), parametrized by a 



dr 



1 



e. 



d<7 

da 3 
i.e., we may write Eq. (|15p as 

dlnn" 2 



da 3 
The solution of Eq. (flH)) . is 



e. 



(16) 



n°(r(a),p(a)) = n"(r(ao),p(ao)) ' -(17) 

where cto corresponds to some initial value of the param- 
eter a. Equation (|17p expresses the fact that the vari- 
ation of the number of NTPs per unit of volume along 
a certain curve (parametrized by a) is directly related 
with the variation of the rest-mass density of the ther- 
mal plasma between the initial and final points of such a 
curve. 

We now turn back to Eq. psp and derive the form 
of the collisions term in the case that the only relevant 
radiative losses are d ue to synchrotron processes . In such 
a case we have (e.g., iRvbicki fc Lightmanlll979f ) 



iarP^Us 



S(p,r) 



(18) 



where ctt is the Thompson cross section, mc is the elec- 
tron rest-mass, J7b = b^/87r is the magnetic energy den- 
sity, and we have assumed that the electrons are ultra- 
relativistic, p ~ ^rric, 7 being the electron Lorentz factor 
(not to be confused with the plasma Lorentz factor). If 
8 = 0, Eq. (|15|) reads in the comoving frame 



i:>lnn 
Dt 



XldntC-l 



(19) 



and, on the other hand, the particle number conservation 
yields 



Dt 

or, equivalently, 

'Dlnn° 



Dt 



syn 



syn 



dB ^ ainn" 
dp dp 



(20) 



Taking into account Eq. p^ . we may plug Eq. (|20p 
into Eq. p5|) to account for the combined effects of syn- 
chrotron losses and adiabatic expansion/compression of 
the fluid 



£)lnn 
Di 







B 



dp 



3 dp 



(21) 



The formal solution of Eq. ([2T|) . can be found following 
the same procedure we used above for the homogeneous 
continuity equation. In this case we interpret 



Din 77.' 
Dt 







P 



Q + B 



9 In 77° 
i9p 



as the derivative of In 77 along the curve 



da 
dp 

da 



B, 



which yields, on the one hand 
dp P 



dT 

and, on the other hand, 



e + s(p,T) 



(22) 



(23) 



n''{T{a),p{a))^n°{T{ao),p{ao)) 



X exp 



da' 



,96(p,r), 



dp 



W) -(24) 



Equation shows the evolution of the particle mo- 
mentum in time, while Eq. (j24p is only a formal solution 
since the exact dependence of p(a), necessary to perform 
the integration, is only known through the differential 
equation ([22|) . The first term on the right hand side of 
Eq. (P^l) accounts for the change of momentum due to 
the adiabatic expansion or compression of the fluid in 
which NTPs are embedded. The time dependence of B 
is fixed by the hydrodynamic properties of the thermal 
fluid and by the comoving magnetic field b, assumed to 
be provided by hydrodynamic simulations and models of 
the b- field (which is not directly simulated) , respectively. 

In order to speed up the numerical evaluation of 
Eqs. ((23)) and (|24|) . we assume that both, the fluid expan- 
sion and the synchrotron losses (or, equivalently, C/b) are 
constant within an small interval of proper time around 
t(cto). Thus, we can write Eq. (22) as 



dp , 

^ = fcaP 

da 



ksP^ 



(25) 



fca and fcg being both constants, such that the following 
relations hold 



p{T{a)) 



f f - (26) 

B{p{a),T{a))^-hp\ (27) 

with Act — a — ao- Equation (|25|) has the following 
analytic solution 



P('7) = Po 



fca + Poks {e^-^" - 1) 



(28) 



where po := p(cto). Upon substitution of the rela- 
tions ([^ and ((Ml) in Eq. we obtain 



n^{T{a),p{a))=n°{T{aQ),p{an)) x 



l + Po 



ka 



r 



(29) 



This equation is approximately valid in the neighbor- 
hood of t(cto) or if the fluid expansion and magnetic field 
energy are both constant in a certain interval Act. In- 
deed, such an assumption is adequate for our purposes, 
since the hydrodynamic evolution is performed numer- 
ically as a succession of finite, but small, time steps. 
Within each hydrodynamic time step the physical vari- 
ables inside of each numerical cell do not change much 
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and, thus, the magnetic field energy and the fluid expan- 
sion are roughly constant. Alternatively, one might not 
assume anything about Q or Ub and solve the system 
of integro-difFcrcntial equations ((22)) and p4l) . However, 
such a procedure is much more computationally demand- 
ing than obtaining the evolution of p and n'^ from, re- 
spectively, Eqs. and Furthermore, since the 
magnetic field is assumed in this work, i.e., not consis- 
tently computed, a numerical solution of the aforemen- 
tioned equations does not yield a true improvement of 
the accuracy. 

For completeness, as in the diffusion approximation 
= 4:7rp^f^^\ we can specify the evolution equation 
for the isotropic part of the distribution function of the 
NTPs 

4 



/(")(T(a),p(a))==/, 



(0) 



1 



r 



,(30) 



where =/W(r(ao),p('^o)). 

Finally, we define the number density of NTPs within 
a certain momentum interval [pa(T(cr)), pt,(T((T))] 

A^(T(a);p„,pb):= ^ "dpn°(T(ao),p(^o)),(31) 
■Jp.iri'yy) 

whose evolution equation can be easily obtained from 
Eqs. dSHD and ^ and reads 

AA(T(a); Pa, pb) = e'^'^^^^AfiTiao); Pa, P^) • (32) 

Eguation ip^ shows that the time evolution of the 
number density of NTPs in a time-evolving momentum 
interval, depends only on the adiabatic changes of the 
NTPs in such momentum interval, but not on the syn- 
chrotron losses (Eq. ([5^ is independent of fcg). 

3.2. Discretization in momentum space 

In the following we normalize p to moC, which allows 
us to express our results in terms of the particle Lorentz 
factor 7 instead of p. In order to make Eqs. and (|29p 
amenable to numerical treatment, we discretizc the mo- 
mentum space in iVb bins, each momentum bin i having a 
lower bound 7^. In the present appl ications we use iVb = 
32 (se c App. IA.4.21). As in. e.g.. iJones. Rvu fc Engel 
(|1999t ). lMiniatil (|200lh . or lJones fc Kand (|2005[ ). we ini- 
tially distribute 7,; logarithmically, i.e. 

(*-l)/(JVb-l) 



7n 



7n 



In 



7min and 7max bciug the minimum and maximum Lorentz 
factors of the considered distribution, respectively. 

On the other hand, the time dimension is also dis- 
cretized in time steps. We call r" the interval of proper 
time elapsed since the beginning of our simulation, and 
denote At = t"+i - r" . 

Our numerical method follows the time evolution of 
NTPs in the momentum space employing a Lagrangian 
approach. We track both the evolution of the A'^b inter- 
face values n° (from Eq. (|29p . where we take t = a and 
CTo = t(o-o) r"). 



n^(r"+^):=n"(r"+\7«(r"+')) =^"(r",7«(r")) x 

n 2 



e'^^-'-(l+7.(r")^(e 
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(33) 



as well as the A^b bin integrated values 



M(7 



7i+l(-r) 



d'j n^{T, 7) . 



(34) 



7i(i') 



The time evolution of the A^b + 1 interface values 7; (r) 
is governed by Eq. 

For the purpose of efficiently computing the syn- 
chrotron emissivity (see §|3]), inside of each Lorentz fac- 
tor bin I, we assume that, at any time, the number of 
NTPs per unit of energy and unit of volume (r, 7) 
< 7 < 7i-i-i (■'')) follows a power-law and, there- 
fore, the whole momentum distribution of NTPs consists 
of a piecewise power-law and. 



nO(T,7) = n^(T)(^ 



1, 



,A^h 



(35) 



where nUr) is the number of particles with 7 = 7^ at 
the proper time r, and qi{T) is the power- law index of 
the distribution at the i-Lorentz factor interval. The val- 
ues of qi{T) are computed numerically in every time step 
plugging Eq. ([55)1 into Eq. and solving iteratively 
the corresponding equation (which also involves knowing 
the interface values -Eq. ([33]) - and justifies why we need 
to follow the evolution of two sets of variables per bin) . 

The approach defined up to here has the advantage 
that, at every time level r", the momentum-space evo- 
lution and the physical space trajectory of the NTPs are 
decoupled during the corresponding time step Ar. The 
hydrodynamic evolution of the thermal plasma provides 
the values of fca and fcg at the beginning of the time step 
(r — r"), and once these values arc known, it is pos- 
sible to compute the momentum distribution of NTPs 
at time r"+^. Thereby, it is possible to perform sepa- 
rately the trajectory integration of the NTPs once, and 
to evolve NTPs in the phase space afterward, as many 
times and with as many initial particle distributions as 
desired (viz., during a post-processing phase). 

3.3. Normalization and initialization of the NTP 
distribution 

Our models are set up such that we initially inject 
through the jet nozzle NTPs with a momentum distri- 
bution function which follows a single power-law, i.e., 
qi = (?i,Vi. Therefore, the initial number and energy 



density in the interval 7n 



< T < 
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7niii: 



2-qi 



(36) 
(37) 



Consistent with our assumptions about the relation be- 
tween the thermal and non-thermal populations we as- 
sume that N = Cj^p/me and U = c„P, where and 
Cj^ are constants, while P and p stand for the pressure 
and rest-mass density of the background fluid, respec- 
tively. Such proportionalities along with Eqs. ()36p and 
([37ll yield (G95) 



7n 



gi - 2 P 1 - (7max/7min)^~'^^ 
- 1 pC? 1 - (7max/7min)^"''i ' 



(38) 
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Fig. 1. — Top panel: Synthetic (total intensity) radio map of the quiescent OP-L-NL model computed at an observational frequency of 
43 GHz. Bottom panel: Same as top panel but for model OP-L-AM, and plotted using the same intensity scale. To compute the model 
OP-L-NL, 32 Lagrangian particles evenly distributed across the jet nozzle have been let to evolve. Since SPEV-NL does not include the 
effect of synchrotron losses on the NTP evolution, the differences between both radio maps are small. In both panels a 10° jet viewing 
angle is assumed. 



and we can use cither Eq. ([55)) or Eq. ([57]) to compute 
nl if the ratio Cj := 7max/7min is fixed. Thus, the initial 
distribution of particles can be determined from pressure 
and rest-mass density at the jet nozzle, simply by speci- 
fying and c„ . 

A key difference between SPEV and AM methods is 
that in SPEV the dimensionless proportionality param- 



case b is ordered, a is the angle the comoving magnetic 
field forms with the line of sight, and vq = 3e|b|/47rmoC. 
In the previous expressions, F is the first synchrotron 
function 



F{x) 



(40) 



, J 1 -c J J- j-i • 4- • • A'c/Q being the modified Bessel function of index 5/3, 

eters c., and c,, are only specified at the jet mjection '^{■^-^ • , ^ i 

^ " ^ J J g^j^^ jl defined as 



nozzle. In the SPEV method, the subsequent time evo- 
lution of the NTP momentum distribution, namely, the 
spectral shape (piecewise power-law) and the limits of 
the distribution 7niin and 7inax as it evolves in the physi- 
cal space is computed according to Eq. (^5]) . AM ignores 
synchrotron loses, which yields a fixed power-law index 
for the whole distribution of Lorentz factors of the NTPs. 
The remaining two parameters needed to specify the dis- 
tribution function, 7min and 7inax are computed from the 
local values of the hydrodynamic variables. On the one 
hand, 7inin follows from Eq. ([55)) and 7max is obtained 
from the fact that, Cy is strictly constant in time if the 
evolution is adiabatic. Also, in contrast to SPEV, it is 
necessary to assume a value of C-y everywhere in the sim- 
ulated region and not only at the injection region. 

4. SYNCHROTRON RADIATION AND SYNTHETIC RADIO 

MAPS 

The synchrotron emissivity. at a time r, of an ensemble 
of NTPs advected by a thermal plasma element, can be 
cast in the following general for m (valid both f or ordered 
and random magnetic fields; see [Mimical I2OOI 



V) 



47rTOoC^ 



E 



d7n°(T,7)5 



where {g{x)^h^^v\}j = (i?(.T), |b|, i/q) if b is randomly 
oriented, or {g(x)^h^^v\}j = (i^(a;), |b| sina, I'o sina) in 



R(x) 



dasin^a f{^— \ . (41) 
V sin a I 



The synchrotron self-absorption process is also in- 
cluded in our algorithm. Thus, we need to compute the 
synchrotron absorption coefficient, at a time r, of an en- 
semble of NTPs advected by a thermal plasma element, 
which can be cast in the following general form 



K(r, v) - 



/•7> + i(t) 

,;=1 



(42) 



d7 



In order to perform the integrals of Eq. ([55]) and ([^^ . 
it is necessary to make some assumption about the inter- 
nal distribution of NTPs within each Lorentz factor bin 
i. As explained in Sect. 13.21 we choose to assume that 
NTPs distribute as power-law (Eq. [35)) inside of each bin. 
This choice agrees with the common assumptions made 
in the literature and is also supported by theoretical ar- 
guments and observations of discrete radio sources (e.g., 
iBicholczvk 1973, cha pt. 6;[Km igl 1 981 and bv numer- 
ical simulations (e.g.. lAchterberg. et al.|[2001f ). Further- 
more, it allows us to build a very efficient and robust 
method for computing the local synchrotron emissivity 
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Fig. 2. — Left and right panels correspond to quiescent models PM-L and PM-H, respectively. The upper panels display the spectral 
energy evolution along the jet axis. Thick black, red and green lines show the evolution of 7max for initial values of equal to 10, 10'^ 
and 10"^, respectively. In the upper panels wo overplot the values of the Lorentz factor corresponding to the maximum emission efficiency 
(solid blue line). The lower and upper 7 where the efficiency drops below 10% of the maximal as a function of distance for an observational 
frequency of 43 GHz are displayed by dotted blue lines above and below the maximum efficiency line, respectively. These blue lines have a 
positive slope since the magnetic field decreases with distance, so that ever larger Lorentz factors are needed to emit efficiently at a given 
frequency. The lower panels show the synchrotron emissivity of each model. The inset in the right panel shows the logarithm of 7max as a 
function of the logarithm of Z for the three models. It can be seen, that for the model PM-H, 7max becomes virtually independent of its 
initial value at a distance larger than 1 i?;, from the jet nozzle. Please, note the difference in the scales of emissivity for the lower panels. 



and the local absorption coefficient, ft consist of tabu- 
lating the functions F{x) and and then tabulating 
integrals over power-law distributions of particles. Pro- 
ceeding in this way is ^ f 00 times faster than computing 
Eqs. (|5^ - (|¥T|) by direct numerical integration. 

The synchrotron coefficients fEqs.[39land l42)) of steady 
jet models result from the time evolution of the La- 
grangian particles injected at the jet nozzle and spatially 
transported along the whole jet (the larger the number 
of Lagrangian particles, the better coverage of the whole 
jet). In our simulations around A'stcady = 32 of such 
Lagrangian particles (i.e., about 4 particles per numeri- 
cal cell at the injection nozzle) are sufficient to properly 
cover a quiescent jet. If the jet is not steady, e.g., be- 
cause a hydrodynamic perturbation is injected, we need 
to follow many more Lagrangian particles. It becomes 
necessary to have particles everywhere the quiescent jet 
is perturbed. For the models in this paper, this means 
to inject new Lagrangian particles through the jet nozzle 
at all time steps after a hydrodynamic perturbation is 
set in. The distribution function of the NTPs injected 
with the perturbation is the same as that of the parti- 
cles injected in the quiescent jet. This is justified since 
the perturbation only changes the bulk Lorentz factor, 
but not the pressure, or the density of the fluid. In the 
simulations where we have injected a hydrodynamic per- 
turbation this implies that we must follow the evolution 
of more than A^stoady x iVtimostcps > 10^ Lagrangian parti- 



cles. This makes our SPEV simulations effectively four- 
dimensional (two spatial, one momentum and a -huge- 
number of Lagrangian particles dimension). Therefore, 
the spatial resolution that we may afford results severely 
limited. 

The synchrotron coefficients depend on the magnetic 
field strength and orientation as well as on the spectral 
energy distribution vP{t^ 7). In our models the magnetic 
field is dynamically negligible, thus we set it up ad hoc. 
We choose that remains a fixed fraction of the particle 
energy density and that the field is randomly oriented. 

Synthetic radio maps are build by integrating the 
transfer equations for synchrotron radiation along rays 
parallel to the line, accounting for the appropriate rel- 
ativistic effects (time dilation, Dopplcr boosting, etc.). 
The technical details relevant for this procedure can be 
found in Appendix [X] 

5. RADIO EMISSION 

The goals of this section are twofold. First, we validate 
the new algorithm comparing the synthetic radio maps 
obtained with SPEV without accounting for synchrotron 
losses with the ones obtained using AM. For this purpose, 
we will employ the SPEV method to evolve NTPs but 
taking fcg = in Eq. ([^5)1 . We will refer to this method of 
evaluating the evolution of NTPs as SPEV-NL. Second, 
we will show the differences induced by accounting for 
synchrotron losses in the evolution of NTPs. 
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Fig. 3. — In the whole figure, black (red) lines correspond to models computed with the SPEV (SPEV-NL) method. The left and 
right panels display properties of the PM-L and OP-L models, respectively. Upper panels: Spectral energy evolution along the jet axis 
for the stationary jet models. The values shown are computed in the jet comoving frame but the distance along the axis is measured 
in the laboratory frame (attached to the jet nozzle). The parameters of the spectral distribution of NTPs are the same as the ones of 
the reference model (§[5]l. In the upper panel the thick lines track the values of 7niin (solid) and 7max (dashed) of SPEV and SPEV-NL 
electron distributions as a function of the distance along the jet axis. Note that the lines corresponding to the values of 7inin fof SPEV and 
SPEV-NL arc almost indistinguishable. The curvature in the line corresponding to 7^ax^> specially in the first 30ii(,, shows the effects of 
synchrotron co olin g of the highest-energy SPEV particles. The blue lines have the same meaning as in Fig. |2] Lower panels: Synchrotron 
emissivity fEg. I39I I at the jet axis is shown as a function of the distance to the jet nozzle. For the parameters chosen, most of the electrons 
of both SPEV and SPEV-NL distributions emit synchrotron radiation efficiently in the whole jet. This makes that both, SPEV-NL and 
SPEV methods display a very similar emissivity dependence with distance along the jet axis (both curves are almost coincident). 
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Fig. 4. — Upper panels: Left and right panels display properties 
(as seen by a distant observer) of the PM-L and OP-L models, 
respectively, computed with the SPEV method. The thin-solid, 
dashed and thick solid lines correspond to the specific intensity 
at frequencies 43 GHz, 22 GHz and 15 GHz respectively. The in- 
tensities are obtained directly from the models without convolving 
the data. For clarity, all the specific intensities are normalized to 
a common value. The dotted line shows the spectral index aia. 
Lower panels: Same as the upper panels, but for the models PM-H 
and OP-H. 



In order to properly compare SPEV-NL and AM re- 
sults we set up the same spectral parameters at the 



jet nozzle for both: qi = 2.2, 



330, 



10^ 



and Pa = 2 X lO^^^gcm"'^. We produce all our im- 
ages for a canonical viewing angle of 10° and assuming 
that i?f, = 0.1 pc. The comoving magnetic field strength 
is ht := Vb2 = 0.02 G (model PM-L-NL) and 0.03 G 
(model 0P-L-NL).4 

For the set of reference parameters we have considered, 
the synthetic radio maps of the quiescent jets produced 
with SPEV-NL yield very small differences with respect 
those computed with AM (Fig. [T|). Indeed, the overall 
agreement between both methods in the predicted qui- 
escent radio maps is remarkably good, particularly, if 
we consider the fact that SPEV is a Lagrangian method 
while AM is Eulerian. 

Looking at the synthetic radio maps of model OP-L- 
NL (Fig. [T]), we observe, as in G97, a regular pattern 
of knots of high emission, associated with the increased 
specific internal energy and rest-mass density of inter- 
nal oblique shocks produced by the initial overpressure 
in this model. The intensity of the knots decreases along 
the jet due to the expansion resulting fr om the gradient in 
external pressure. Some authors (^e.g.. lDalv fc Marscheil 



5.1. Calibration of the method 



^ With such values of b(, the magnetic field is dynamically neg- 
ligible. 
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Fig. 5. — The number of particles per unit of energy and unit 
of volume n^{-f) is displayed with a solid black line as a function 
of 7 for the P M-L model. We omit the temporal dependence of 
nP{T, 7) in Eq. Il35|l because we are considering quiescent jet mod- 
els. The solid blue and red linos show the synchrotron function 
Rv{x) (Eg. I4ll l at frequencies 15 GHz and 43 GHz, respectively, 
while the products n'^('y)Ri5{x) and ".'^(7)^43(2:) arc displayed 
with dashed blue and red lines, respectively. The later pro duct s 
are precisely the integrand of the synchrotron emissivity (Eq. I39I I. 
Lower panel: Corresponds to the typical conditions one encounters 
downstream the jet (7min ~ 135, 7max ~ 6 X 10^, b ~ 0.002 G). 
Upper panel: Corresponds to the conditions found close to the in- 
jection nozzle (7min = 330, 7max = 3.3 X 10^, h ~ 0.02 G). 



an 
a 



:r\ PM-L 

.1.1.1. i^i*?- 


OP- 1, 

^ V/ \ 
/ - ^ \ 

1,1,1, 1 "i^r'^Rv 












/ ^ ^ ^\ /^ 




— --A^ 


5 10 15 20 25 


5 10 15 20 25 
^obsIRJ 



-0,2 



■0,4' 



■0,6 



■2 

■2,5 



Fig. 6. — Same as Fig. |4] but in this case the specific intensi- 
ties are obtained by convolving the data with a circular Gaussian 
beam whose radius at FWHM are 2.25R(,, 4.40if,5 and %A^R^, at 
frequencies 43 GHz, 22 GHz and 15 GHz, respectively. With this 
convolution we degrade the resolution of our data to limits compa- 
rable with VLBI observations. 



Il988t G95; G97; iMarscher. et al.|[2008l ) propose that the 
VLBI cores may actually correspond to the first of such 
recollimation shocks. Since, for the parameters we use, 
the source absorption for frequencies above 1 GHz is neg- 
ligible, the jet core reflects the ad hoc jet inlet in the 
PM-L-NL model, while we shall associate it with the first 
recollimation shock for model OP-L-NL. The rest of the 
knots are standing features in the radio maps for which. 
there exists robust observational confirmation (|G6mea 
[200l[200l . 




10 10 10 
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Fig. 7. — The different lines in the plot show the spectral energy 
distribution of selected points along the jet axis of model PM with 
two diff'erent magnetic fields. All the models have been computed 
using the SPEV method with the reference parameters of §[5] Solid 
lines correspond to the PM-L model and dashed lines to the PM-S 
model. The distance to the nozzle (in Kh) to which each spectrum 
corresponds is provided in the legend. Synchrotron self-absorption 
is dominant at frequencies below few hundred MHz. 



Since the synchrotron losses affect more the higher en- 
ergy part of the distribution of NTPs than the lower 
one, we have also validated our code by considering 
the dependence of the results with the limit 7niax and 
checked them ag ainst the theoretical expectations (e.g., 
lPacholcz"v3ll970[ ). For this we reduce the value of 7niax 
keeping all other parameters fixed and equal to those of 
the PM-L model. Since the value of 7max is set by the 
ratio C^, in order to study the dependence of the re- 
sults with 7max, we have computed a set of models com- 
bining three different values of = {10"^, 10^, 10} and 
bft = 0.02 G. Additionally, to highlight the effect of the 
radiative losses, we have performed the same simulations 
(varying C^,) for a larger value of the beam magnetic 
field, equal to that of the model PM-H. 

For model PM-L (Fig. [21 left panel), radiative losses are 
negligible, and the reduction in (i.e., in 7niax), does 
not change appreciably the radio maps at radio observing 
frequencies. Indeed, except the model with the lowest 
(corresponding to 7max = 2200) beyond 160 all 
the models stay above the 100% efficient radiation limit 
along the whole jet. 

The models with larger magnetic field b^ = 0.2 G 
(Fig. [21 right panel), undergo a much faster evolution. 
The emissivity along the jet axis drops very quickly and 
at z = 150i?b, it is five orders of magnitude smaller 
than for the PM-L model. After a very short distance 
(~ 1 i?b), synchrotron losses bring 7max of all three mod- 
els to a common value which is independent of the initial 
one (note that the variation of 7max with distance is in- 
distinguishable for the three models except in the zoom 
displayed in the inset of the top right panel of Fig. [2]). 
The reason for this degenerate evolution re sides in the 
relati vely large magnetic field strength (see iPacholczvkl 
[19701 Eq. 6.20). Thus, our method is able to reproduce 
the common evolution of models with different values of 
7max and a relatively large magnetic field. 

5.2. On the relevance of synchrotron losses 
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Fig. 8. — Same as Fig.[T] but in this case for the OP-H model. 
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Fig. 9. — Same as Fig.[3]but for the PM-H and OP-H models. In this case, the synchrotron losses in SPEV are so important that, 7max^ 
leaves the efficient synchrotron radiation regime. The point where this happens depends on the model. The crossing of the 7m^x^ with 
the line denoting a 10% synchrotron efficiency limit occurs when the particles reach z ~ 50ii(, for model PM-H and much earlier z ~ 24/?;, 
in case of model OP-H. After this line crossings the whole distribution radiates very little at the considered observational frequency. This 
produces a substantial dimming of SPEV models at large z. The difference in the emissivity as a function of the distance z is larger for 
the PM-H model than for the OP-H model because of the re-compressions that the SPEV-NTP experiences at shocks. 



Having verified that our method (SPEV-NL) compares 
adequately to the AM, we now turn to the specific role 
that synchrotron losses play in the evolution of NTPs. 
For that, we compare in Fig. [3] the spectral properties of 
NTPs in quiescent jet models using both SPEV-NL and 
SPEV methods. It is obvious that the highest energy 



particles of the distribution cool down rather quickly (see 
the fast decay of the dashed black curves in the upper 
panels of Fig. ^ even for the small value of b^ consid- 
ered here. Most of the spectral evolution triggered by 
synchrotron cooling at high values of 7 happens in the 
first 25i?b — 50i?f,. After that location, the ratio is 
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much smaller than at the injection nozzle (C^ < 50), and 
the evolution of the NTP population is dominated by 
the adiabatic cooling/compression downstream the jet. 
In contrast, the upper limit of the SPEV-NL distribu- 
tion (dashed red curves in the upper panels of Fig. [3]) 
only changes by a factor of 2 along the whole jet length. 
Theoretically, it is well understood that it is possible to 
undergo a substantial spectral evolution (triggered by 
synchrotron losses) and, simultaneously, not to have any 
manifestation of such evolution at radio frequencies (e.g., 
lPacholczvklll970f ). The substantial decrease of 7max trig- 
gered by the radiative losses, does not affect much the 
value of the integral that has to be performed over 7 in 
order to compute the emissivity in Eq. (|39p . since most 
of the emitted power at radio-frequencies happens rela- 
tively close to 7min, whcrc synchrotron losses are negligi- 
ble. Certainly, at higher observing frequencies this is not 
the case, and the emissivity substantially drops because 
of the fact that both, the synchrotron losses (Eq. [T8|) 
and the frequency at which the spectral maximum emis- 
sion is reached depend on the square of the non-thermal 
electron energy (and on the magnetic field strength). 

We define the spectral index between two radio fre- 
quencies as 

_ logiS,/Sj) 
"'^ log(zy,/:/j) ' 



(43) 



where Si and Sj are the flux densities at the frequencies 
Vi and Vj , respectively. Since we compute synthetic radio 
maps at three different radio frequencies {vi = 15 GHz, 
1/2 = 22 GHz, and 1/3 — 43 GHz), we may define three 
different spectral indices. For convenience, in the follow- 
ing, we consider the spectral index ai3 between 15 GHz 
and 43 GHz. Furthermore, we may compute aia for both 
convolved or unconvolved&ux densities. The unconvolved 
flux density is directly obtained from the simulations and 
has an extremely good spatial resolution, viz. the uncon- 
volved radio images have a resolution comparable to that 
of the hydrodynamic data. The convolved flux densi- 
ties result from the convolution with a circular Gaussian 
beam of the unconvolved data. The FWHM of the Gaus- 
sian beam is proportional to the observing wavelength. 
This convolution is necessary to degrade the resolution 
of our models down to limits comparable with typical 
VLBI observing resolution. We note that in order to 
compute the spectral index for convolved flux densities, 
we have to employ the same FWHM convolution kernel 
for the data at the two frequencies under consideration. 
Thus, to compute ai3 for convolved data, we employ the 
same Gaussian beam with a FWHM 6.45 Rb for both flux 
densities at 15 GHz and 43 GHz. 

Our models are computed for an electron spectral in- 
dex q = qi = 2.2. We verify that, at large distances to 
the jet nozzle, unconvolved models (Fig. [4] upper panels) 
tend to reach the expected value a = {1 — q) /2 = —0.6 
for an optically thin source. This asymptotic value is 
reached smoothly in case of the PM-L and PM-H models 
and it is modulated by the presence of inhomogeneities 
(recollimation shocks) in the beam of models OP-L and 
OP-H. 

Close to the jet nozzle, our unconvolved models display 
flat or even inverted (ai3 > 0) spectra (Fig.|4]), in spite of 
the fact that the jets are optically thin throughout their 
whole volume. The occurrence of flat or inverted spec- 



trum depends on the magnetic field strength and differs 
for OP and PM models. As shown in Fig. [H the PM-L 
model shows an inverted spectrum for z < 2.5Rb, while 
the OP-L model displays a pattern of alternated inverted 
and normal (ai3 < 0) spectra for z < 12.5i?f,. The spec- 
tral inversion in the OP-L model happens where standing 
features (associated to recollimation shocks in the beam) 
are seen in the jet. 

If synchrotron losses are not included, the spectral be- 
havior of models PM-L and OP-L remains unchanged, 
because in such a case loses are negligible. However, if 
for the models PM-H and OP-H the losses are not ac- 
counted for (which is, obviously, a wrong assumption), 
the jet displays an inverted spectrum up to distances 
z - 30i?fc. 

The behavior of the spectral index exhibited by our 
models close to the jet nozzle contrasts with the theoret- 
ical expectations for an inhomogeneous, optically thin 
jet with a negative electron spectral index, for which the 
jet in homogeneity is predicte d to steepen the spectrum 
(e.g.. lMarscIierlll980l : [Konig]|ll98lD . To explain this dis- 
crepancy we argue that the analytic predictions are based 
on the assumption that the limits of the energy distribu- 
tion of the NTPs safely yield that the contribution of 
the synchrotron functions (Eqs. |40] and l4T]l to the syn- 
chrotron coefficients (Eqs. [39l and [42|) is proportional to 
some power of the frequency and of the NTP's Lorentz 
factor. This situation does not happen if the lower limit 
of the distribution nP{j), 7min, is (roughly) larger than 
the value 7,^, at which the synchrotron function i?(a;iow) 
(Eq. |4T|) reaches its maximum, where xiow = i^iow / i^ol'^ 1 
and flow is the smallest observing frequency in the co- 
moving frame. Since the function R{x) has a maximum 
for X ~ 0.28, one finds that the condition to have an in- 
verted spectrum is 7min > 7m - 1-9 ^ (yi^^ Ivq^I'^V^^I'^ , 
where T) := l/r(l — ucos0) is the Doppler factor. Since, 
in our case, v\ov, = 15 GHz, we may also write 



7min > 113 



\l/2 f h_ 

15 GHz/ VTg 



-1/2 



V 



-1/2 



(44) 



Figure [5] shows how this boundary effect substantially 
modifies the emissivity at 15 GHz and 43 GHz for the 
model PM-L. At the injection nozzle (Fig.[5]upper panel) 
the lower limit of the integral in Eq. ([55]) is set by 7niin 
and not by the lower limit of Ryix). However, down- 
stream the jet (Fig. [5] lower panel) the situation reverses 
and the fast decay of Ru{x) for 7 < 300 sets the lower 
limit of the emissivity integral. Thus, close to the nozzle, 
the value of the area below the nP{'j)R43{x) curve, which 
is proportional to the emissivity at 43 GHz, is larger than 
that below the curve n'^{j)Ri5{x). Hence, there is an 
emissivity excess at 43 GHz compared to that at 15 GHz. 
As a result, the 0:13 becomes positive close to the jet noz- 
zle. Far away from the nozzle the emissivity at 15 GHz 
almost doubles that at 43 GHz, explaining why values of 
q;i3 < are reached asymptotically. 

The convolved models display some traces of the be- 
havior shown for the uncovolved ones. For example, OP 
models display a flat or inverted spectrum very close to 
the jet nozzle (Fig. [5] right panels). This is not the case 
for PM-L model (Fig.lHl). Since the resolution of the con- 
volved data is much poorer than that of the unconvolved 
one, ai3 exhibits a quasi monotonically decreasing profile 
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from the jet nozzle (where —0.1 < aia < 0). The coarse 
resolution of the convolved data also blurs any signature 
in the spectral index associated to the existence of cross 
shocks in the beam of OP models. Furthermore, the de- 
cay with distance of the spectral index is shallower for the 
convolved flux data than for the unconvolved one. Hence, 
the theoretical value ais = —0.6, which is expected for 
an optically thin synchrotron source, is reached nowhere 
in the jet models PM-L and OP-L (Fig. 

As expected, at frequencies below a few hundred MHz, 
the jet is strongly self-absorbed everywhere (Fig. [7]). 
Close to the jet nozzle, there is not a clear turnover fre- 
quency between the self-absorbed part of the spectrum 
and the optically thin one. Instead, we observe a smooth 
transition between both regimes. Far from the nozzle, 
the self -absorption turnover is much more peaked. It is 
known (jTsang fc Kirkll2b07f ) that in contrast with a dis- 
tribution of NTP that follows a power-law extending to 
7min ^ 1) if the power-law is restricted to a relatively 
large, but not unrealistic 7min, or if the electron dis- 
tribution was monoenergetic, the intensity can be flat 
over nearly two decades in frequency (which implies that 
the energy flux grows linearly over the same frequency 
range). Our PM-L models have 7,„i,i — 330 at the in- 
jection nozzle and reduce it to 7min — 200 at z = 200i?6 
because of the adiabatic expansion of the jet (Fig. [3] up- 
per left). As we have argued in §[S]close to the jet nozzle, 
7min > 7m j which means that 7min is sufficiently large to 
be in the range where a sm ooth turnover t ransit ion is 
expected, in agreement with iTsang fc Kir3 (|2007[ ). Far 
away from the nozzle, since 7min decreases, we recover 
the more standard situation in which an obvious turnover 
frequency can be identified. 

Provided that close to the nozzle our PM (also OP) 
models are weakly self-absorbed (at 15 GHz, the solid 
black line in Fig. [7] has not reached the power-law regime 
yet), one may question whether the spectral inversion we 
have found is not also the result of opacity effects. We 
have dismissed such a possibility by running models with 
the SPEV method including no absorption. 

5.2.1. Dependence with the magnetic field strength 

In order to study the effect of intense synchrotron 
losses we consider models PM-H and OP-H (Figs. [8] and 
[9]). Very close to the injection nozzle {Z ^ 50 Rb) the 
line denoting the evolution of 7max^ crosses the line cor- 
responding to a 10% synchrotron efficiency limit (lower 
blue thick line; Fig.lH]) and most of the synchrotron emis- 
sivity falls outside of the observational frequency. Be- 
cause of a stronger magnetic field than in models PM-L 
and OP-L, more energy is lost close to the jet nozzle 
than far from it and, thus, SPEV-radio-maps look much 
shorter than SPEV-NL-radio-maps (Fig. [8]). An alterna- 
tive way to see such an effect is through the rapid decay of 
7max^ in the first 10i?b in Fig. [21 right panel. Afterwards, 
the adiabatic changes dominate the NTP evolution. The 
initial period of fast evolution is even shorter if a larger 
magnetic field were to be considered. 

The intensity contrast between shocked and unshocked 
jet regions of model OP-H (Fig. [9]) is larger than that 
of model OP-H-NL. Indeed, the OP-H model appears 
as a discontinuous jet (Fig. [H]) because of the slightly 
larger intensity increase than in the OP-H-NL model 



when the NTP distribution passes through cross-shocks 
and the much more pronounced intensity decrease at rar- 
efactions. We note that, although the adiabatic evolu- 
tion is followed with the same algorithm in SPEV and 
SPEV-NL, the radiative losses change substantially the 
NTP distribution that it is injected through the nozzle 
after very short distances. The consequence being that 
the NTP distribution (r, 7) that faces shocks and rar- 
efactions downstream the nozzle is rather different when 
using SPEV or SPEV-NL method and, therefore, the rel- 
ative intensity of shocked and unshocked regions is also 
different depending on whether synchrotron losses are in- 
cluded or not in the calculation. 

The outlined differences between models OP-H and 
OP-H-NL (with shocks in the beam), have to be inter- 
preted with caution since none of the methods accounts 
for the injection of high-energy particles at shocks. But 
independent of this, we expect that if the magnetic field 
is sufficiently large, the SPEV method will yield a rather 
fast evolution of such particles and, thereby, a faster de- 
cay of the intensity downstream the shock. 

The most relevant difference between the upper and 
lower panels of Fig. [8l is the increased brightness of the 
jet close to the injection nozzle and the steeper fading of 
the jet when energy losses are included. This fact poses 
the paradox that the method that accounts for radiative 
losses (SPEV) yields brighter standing features close to 
the injection nozzle (far from the nozzle the situation re- 
verses and the SPEV-NL model is brighter than SPEV 
one). In order to disentangle this apparent contradic- 
tion, we shall consider that the plasma is compressed at 
standing shocks, which yields a growth of the magnetic 
field energy density (proportional to the pressure in our 
case), and triggers a faster cooling of the high-energy 
particles. Since the SPEV method conserves the number 
density of NTPs fEq. [36]) . due to the synchrotron losses, 
high-energy particles reduce their energy and accumulate 
into an interval of Lorentz factor which is smaller than 
in the case of SPEV-NL models. As in such reduced 
Lorentz factor interval NTPs radiate more efficiently at 
the considered radio-frequencies, the emissivity of SPEV 
models at strong compressions (like, e.g., the considered 
cross shocks) becomes larger than that corresponding to 
models which do not include synchrotron losses. It is im- 
portant to note that this situation happens in our mod- 
els rather close to the jet nozzle. The reason being that 
after the NTPs have suffered a substantial synchrotron 
cooling, the evolution of the NTP distribution is domi- 
nated by the adiabatic terms of Eq.[25l In such a regime, 
reached by our models at a certain distance from the jet 
nozzle, the evolution of SPEV-NL and SPEV models is 
qualitatively similar. Considering the different qualita- 
tive evolution of the NTP distribution close to the noz- 
zle and far from it, we refer to such epochs as losses- 
dominated and adiabatic regimes, respectively. These 
terms agree with the commonly used i n the literature to 
refer to similar phenomenologies (e.g.. iMarscher fc Geail 
fT985h . 

For PM-H and OP-H models, the spectral behavior is 
dominated by the change of slope of the NTP Lorentz 
factor distribution beyond the synchrotron cooling break 
at 7 = 7br. Theoretically, an optically thin inhomoge- 
neous jet shall display a spectral index a = {q -\- l)/2, if 
the radiation in the observational band is dominated by 
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the electrons with Lorentz factors 7 > 7br, or a ~ —2.7 if 
the emission is do minated by e lectrons with Lorentz fac- 
tors close to 7max (jK5nig]||1981f )^ . Figure[3] (lower panels) 
shows that asymptotically (viz., at large z) imconvolved 
models reach values ais < — 2.5, implying that the high- 
est energy electrons with 7 ~ 7max are the most effi- 
ciently radiating at the considered observing frequencies. 
The value of 7max differs significantly when synchrotron 
loses are not included. This fact explains the inversion 
of the spectrum along the whole jet if synchrotron losses 
were not included (PM-H-NL and OP-H-NL models). 
Thereby, synchrotron losses tend to produce a "normal" 
spectrum (a^ < 0) if the magnetic field is large. 

Regarding the convolved data, we note that models 
with a higher magnetic field display the same qualitative 
phenomenology discussed in § [51 In this case, the theo- 
retical value a ~ —2.7 is not reached neither by the PM- 
H {oF,f-Z = -1.4) nor by the OP-H (aOP-" = -1.1) 
model (Fig. [6] lower panels). 

6. INFRARED TO X-RAYS EMISSION 

We have computed the spectral properties of some of 
our quiescent jet models above radio frequencies. We 
note, that we have not included any particle accelera- 
tion process at shocks in the SPEV method, thus, the 
spectrum beyond infrared frequencies has to be taken 
carefully. If any shock acceleration mechanism were in- 
cluded, a larger contribution of the shocked regions will 
be present. In addition, the inverse Compton process, 
may shape the emission at such high energies, and such 
cooling process is presently not included in SPEV. 

The results for models PM-S and PM-L (Tab.[T]), which 
have no or extremely weak shocks are displayed in Fig.[7l 
where we show the spectral energy distribution at se- 
lected distances from the nozzle for points located along 
the jet axis. The small magnetic field of model PM-S 
(Fig. [7] dashed lines) minimizes the energy losses, but 
also the observed flux in the optical or X-ray band, ren- 
dering observable at such wavelengths the hydrodynamic 
jet models considered here (if the jet is sufficiently close). 
In the case of model PM-L, right at the nozzle (z = in 
Fig. [7]), the energy flux cut-off is located at ~ 10^* Hz. 
This means that, we could observe the jet core in the soft 
X-ray band, if the source was sufficiently close. How- 
ever, the core size at s uch frequencies is very small (as 
it is expected; see e.g.. iMarscher fc Gearlll985h . This is 
reproduced in our models since at such a short distance 
as z = hRhi the jet can barely be observed in the Near 
Ultraviolet or, perhaps in the optical band (Fig. [71 red 
solid line) , but there is virtually no flux in the X-ray band 
because of the fast NTP cooling for the considered mag- 
netic field energy density at the jet nozzle. In the Near 
Infrared range, the jet could perhaps be observable up to 
distances of 10i?6 — 15i?b. A larger magnetic field drives 
a faster cooling, rendering undetectable the jet even at 
infrared frequencies. This phenomenology has been in- 
voked to explain the relative paucity of optical jets with 
respect to radio jets. However, there are a number of 

^ We obtain t his value from the expression cxs-i = {m + 2 — n)/m 
of IKonigil II1981I ) with m = 1.15 and n = 0. The values of m and n 
are computed from the decay with the distance to the jet nozzle of 
the magnetic field strength |b| oc z~"^ and of the number density 
of NTPs per unit of energy oc z^", respectively. 



authors which claim that a large proportion of jets gen- 
erate significa nt levels of both optic al and, even. X-ray 
emission (e.g., iPerlman. et al.l [20061 ). Our results shall 
not be taken in support of any of the two thesis since 
energy losses depend also on the magnetic field strength 
(Eq. lisp, which we fix ad hoc. 

7. EVOLUTION OF A SUPERLUMINAL COMPONENT 

In this section we discuss the time dependent observed 
emission once a hydrodynamic perturbation is injected 
into the jet (see Sec. 12. 1|) . Following the convention of 
G97, we call components to local increases of the spe- 
cific intensity in a radio map, while we use perturbation 
to denominate the variation of the hydrodynamic condi- 
tions injected through the jet nozzle. In order to mag- 
nify the effect of synchrotron losses in our models, we 
discuss models PM-H and OP-H in Sec. 17.11 and we also 
look for the differences between the PM and OP mod- 
els.^ While the standing shocks of the beam of model 
OP-H are very weak, the shocks developed by the hydro- 
dynamic perturbation are rather strong. Since we have 
not included in our method the acceleration of NTPs at 
relativistic shocks, computing the synchrotron emission 
at frequencies above radio may yield inconsistent results. 
Therefore, we only analyze the spectral properties of the 
emission in radio bands. Finally we show spacetime plots 
of hydrodynamic and emission properties along the jet 
axis in Sec. 17.21 

7.1. On the relevance of synchrotron losses 

The magnetic field energy density is set ad hoc in our 
models (Sect.H]), and we can change it freely if the result- 
ing magnetic field does not become dynamically relevant. 
For the sake of a better illustration of the effect of the 
synchrotron losses on the morphologies displayed in the 
radio maps, we have computed models PM-H and OP- 
H (Fig. [TO]), and PM-H-NL and OP-H-NL (Fig. ]T^. A 
noticeable general characteristic of SPEV-NL models is 
that all the features identifiable in the radio maps are 
more elongated (along the jet axial direction) than in 
the case that synchrotron losses are included. The rea- 
son is that without synchrotron losses, the beam of the 
jet is brighter at longer distances. Thus, in the uncon- 
volved data, the parts located downstream the jet weight 
more in the convolution beam than in the case where 
synchrotron losses are included, biasing the isocontours 
of flux density along the axial, downstream jet direc- 
tion. For the same reason, the models which include syn- 
chrotron losses display a more knotty morphology than 
those which do not include them, both in the uncon- 
volved and in the convolved data. This feature is more 
important in case of OP-H and OP-H-NL models (com- 
pare, e.g., panels two, three and six -from top- of Figs. fTUI 
and [11]) than in case of PM-H and PM-H-NL models. 

The main component undergoes losses-dominated 
(first) and adiabatic (later) regimes as quiescent jet mod- 
els do. In the losses-dominated regime (upper two panels 
of Fig. [T0|) . SPEV models exhibit a brighter component 
than SPEV-NL models. Later, in the adiabatic regime, 
SPEV models display a dimmer component than SPEV- 

^ In the online material we provide a movie ("PMOP- 
fiduc.mpg") where the evolution of the total intensity at 43 GHz is 
displayed for models PM-L and OP-L. 
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Fig. 10. — Snapshots of the emission at 43 GHz due to component evolution computed with the SPEV method, for the PM-H (left) and 
OP-H (right) models. From top to bottom panels show the observed emission 0.02, 0.39, 0.75, 1.12, 1.94 and 4.58 years after the component 
appears. The same gray scale has been used for all snapshots. The superimposed contours have been obtained by convolving the image 
with a circular Gaussian beam whose radius at FWHM is 2.25^6. The contour levels are 0.005, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64 and 
0.9 of the maximum of the convolved emission. The horizontal length scale is expressed in units of Ri, = 0.1 pc, while the vertical length 
scale has been compressed and spans only lOi?,;,. The main component has moved out of the right boundary in the lower two panels. 



NL ones. As we argued in Sect. l5.2?Tl the conservation of 
the NTPs number density explains this phenomenology. 

The main component clearly splits into two parts when 
synchrotron losses are included in model OP-H (Fig. [TUl 
panels 2 and 3 from top; see also the movie "PMOP- 
highB.mpg" in the online material). The component 
splitting is not so apparent in model OP-H-NL. although 
it also takes place farther away from the nozzle than in 
the model including losses (Fig.[TTl third panel from top). 
The splitting of the main component happens during the 
losses-dominated regime and the rear part of the com- 
ponent is brighter than the forward one if losses are in- 
cluded, otherwise, the forward part of the component is 
brighter than the rear one. However, the fact that the 
component is seen as a double peaked structure is not 
the direct result of the splitting of the hydrodynamic 
perturbation in two parts 12. ip . because the projected 
separation of these two hydrodynamic features is smaller 
than the convolution beam, even at 43 GHz. Instead, this 
results from the interaction of the hydrodynamic per- 
turbation with the cross shocks in the beam of model 
OP. Because of the small viewing angle, the increased 
emission triggered in the component when it crosses over 
a recollimation shock is seen by the observer to arrive 
simultaneously with the radiation emitted when the hy- 
drodynamic perturbation was crossing over the preceding 
(upstream) cross shock. 

Figure [T^ shows the evolution of the component at 



15 Ghz (left panels) and 22 GHz (right panels) for the 
PM-H model. The convolution beam depends linearly 
on the wavelength of observation, thereby, it is larger at 
smaller frequencies. Except for the obvious disparity of 
resolutions the evolution of the main component along 
the pressure matched jet at 15 GHz, 22 GHz and 43 GHz 
does not display large differences. The main component 
appears as a moving bright spot at all three frequencies 
(upper three panels of Fig. [10] left and Fig. [T2|). 

We have also checked that the profile outlined above 
does not depend on including synchrotron losses either. 
However, the smaller the magnetic field, the larger the 
increase in the spectral index behind the intensity max- 
ima associated to the main component (i.e., associated 
with the rarefaction trailing the main hydrodynamic per- 
turbation). The time evolution of the prototype spec- 
tral profile of a hydrodynamic perturbation injected at 
the nozzle is characterized by a substantial steepening 
of the spectrum behind the intensity maxima (Figs. 113b 
and [Hb.d) compared to the quiescent jet model. This 
behavior of the spectral index has also been found in 
previous theoretical papers, and it is attributed to the 
fact that the NTP distribution evolves on timescales 
smaller than the ligh t crossing time of the source (e.g.. 
IChiaberee fc Ghiselh ni 1999). 

Comparing Figs. [T3b and [Tib , it is remarkable that 
trailing components pop up precisely to the left (i.e., be- 
hind) of the local relative spectral index minimum (at 
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Fig. 11. — Same as Fig. 1101 but without including synchrotron losses (SPEV-NL method). 
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Fig. 12. — Same as Fig. 1101 but only showing the PM-H model observed at 15 GHz (left) and 22 GHz (right). The images at each 
frequency use an intensity gray scale separately normalized to the maximum at the corresponding frequency. The vertical scale spans 20/?;,, 
i.e., there is a factor of two difference between the vertical scales shown in Fig. 1101 and in this figure. Since the FWHM of the convolution 
beam depends linearly on the wavelength of observation, at 15 GHz and 22 GHz the FWHM are 6.45-Rb and AAORi,, respectively. 
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^obs ^^b-^ ^obs ^^b-^ 

Fig. 13. — In each panel lines correspond to the differences in 
intensity or spectral index between model PM-H with and without 
a hydrodynamic perturbation. Each panel correspond s to a differ- 
ent observer's time (the times are the same as in Figs. [T0l and ll2(l . 
The thick-solid, dashed-dotted and dashed lines represent the nor- 
malized difference (/J'(Z) -/(Z))/ max^ at 15 GHz, 
22 GHz and 43 GHz, respectively. I{Z) and IP{Z) are the intensi- 
ties averaged over cross-sections of the jet at each distance Z from 
the injection nozzle for the quiescent and perturbed models, re- 
spectively. The maximum in the denominator extends for all Z 
along the jet axis. The thin solid line shows the difference in the 
spectral index between the PM-H model with the hydrodynamic 
perturbation and the corresponding quiescent model. Precisely, 
the line shows the function 5 X {ai3{Z) — a^g(Z))/ max^ |cn3(.Z)|, 
where a^^i^) ^-nd ai3{Z) correspond to the cross-sectional average 
of the spectral index (Eg. 143 II of the jet with the injected hydrody- 
namic perturbation and to the quiescent jet, respectively. 
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Fig. 14.— Same as FiglT3]but for the model PM-L. 
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Fig. 15.— Same as FigUHlbut for the model OP-L. 



background jet of the trailing components identifiable at 
43 GHz, depends on the strength of the initial magnetic 
field, in spite of the fact that in our models the mag- 
netic field is dynamically negligible.^ At higher magnetic 
field strength the intensity of the trailing components is 
lowered and, some of them are hardly visible (e.g., the 
leading trailing at ~ 25i?h is evident in Fig. [T4t . while 
it is difficult to identify in Fig. [13?) . Thereby, the ob- 
servational imprint of trailing components is frequency 
dependent. 

The evolution of the perturbation in model OP- 
L displays a slightly different profile at 43 GHz than 
in model PM-L. The main component splits into two 
sub-components at the highest observing frequency 
(Fig. [Hb). At 15 GHz and 22 GHz, the profile of the 
perturbation is qualitatively the same as for the PM-L 
model. The spectral index displays a behavior very sim- 
ilar to that of the PM-L model. However, the evolution 
after the passage of the main component in model OP-L 
(Fig. [15H - [15f ) is different from that of model PM-L. 
The number of bright spots popping up in the wake of 
the main perturbation is smaller and they are brighter 
(in relation to the quiescent jet) in the OP-L model than 
in the PM-L one. Identifying these features as trailing 
components (Sect. [7T2|l . we realize that they do not only 
appear at 43 GHz, but also at 22 GHz, and one may guess 
them even at 15 GHz. 

7.2. Spacetime analysis 

In order to relate the hydrodynamic evolution with 
the features observed in the synthetic radio maps, we 
have built up several space-time diagrams of the evo- 
lution of the component as seen by a distant observer. 
In Fig. [16] we plot the difference in intensity, averaged 
over the beam cross-section, between the perturbed and 



Z - 18i?h in Fig. [HJ; and Z ~ 20i?t in Fig. \T^) that 
follows the local relative maximum of the spectral in- 
dex reached in the wake of the main perturbation. Fur- 
thermore, we notice that the intensity relative to the 



^ According to IMimica. Alov fc Miilleii l|20071 '). the boundary 
separating magnetic fields dynamically relevant from those in which 
the magnetic field is dynamically negligible is around ~ 0.03P. 
In our case, even for the model with the largest comoving magnetic 
field, we have C/b = O.OIP. 
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Fig. 16. — Spacetime plot of the difference in the intensity at 43 GHz, averaged over the beam cross-section, between the perturbed and 
quiescent SPEV emission for PM-L and OP-L (upper left and lower left panels, respectively) and PM-H and OP-H (upper right and lower 
right panels, respectively) unconvolved models. The slope of the dashed line corresponds to an apparent velocity equal to the speed of light. 
The solid black dots correspond to the world-lines of a number of bright features observed in the convolved 43 GHz-radio images resulting 
from the difference between the hydrodynamic models with and without an injected perturbation. Among these features, there are trailing 
components (particularly in the PM models) and also standing recoUimation shocks (characteristic of the OP models). The peaks in the 
color plot (yellow-white shades) do not always match the distribution of black dots due to the difference in the resolution of the convolved 
and unconvolved data. The color scale is linear and common for each column of panels. It displays the difference of averaged intensities in 
arbitrary units. 



quiescent models. This difference accounts for the net 
effects that the passage of the hydrodynamic perturba- 
tion triggers on the quiescent jet. The trajectory of the 
main component is seen as a bright (yellow) region close 
to the top of each plot. Its superluminal motion is appar- 
ent when the slope of the trajectory is compared to that 
of the dashed line, which denotes the slope correspond- 
ing to the speed of light. Below the main component, the 
dark (blue) region is associated to the reduced intensity 
that the rarefaction trailing the hydrodynamic perturba- 
tion leaves. 

As in G97, while in models PM-L and PM-H the main 
component and the reduced intensity region trailing it 
are continuous in the space-time diagrams (Fig. [16] up- 
per panels), in OP-L and OP-H models they flash in- 
termittently as they cross over standing cross shocks of 
the beam (larger intensity ~ Fig. [TBI lower left panel) and 
then expand in the rarefactions that follow such stand- 
ing shocks (smaller intensity). The interaction of the 
perturbation with the standing shocks of the quiescent 
OP model results in a displacement of the position of the 
shocks also noticed in G97. The temporarily dragging of 
standing components, is clearly visible in the lower left 
panel of Fig. [THj The second (from the left) of the well 
identified bright spots, oscillates with an amplitude of 
~ 1.4i?f, in ^ 10 months. The trend being to increase 
both the oscillation period and the amplitude with the 
distance to the jet nozzle. 

Besides the m ain component, we observe several trail- 
ing components ([Agudo et al.ll200"l[ ). identified in Fig. [TBI 
by "threads" with an intensity larger than in the qui- 
escent model, which emerge from the wake of the main 
component. In Fig.[Tn]we also overplot (black dots) the 
world- lines of a number of bright features observed in the 
convolved 43 GHz-radio images resulting from the differ- 



ence between the hydrodynamic models with and with- 
out an injected perturbation. These world-lines show 
only those local intensity maxima which could be unam- 
biguously tracked in convolved radio maps. Except for 
the bright features closer to the jet nozzle, the world-lines 
match fairly well the unconvolved trails of high intensity. 
The latest three trailing components of Fig. [16] (upper 
left panel) do actually recede® in the convolved 43 GHz 
maps as much as 0.5Rb for 1 to 4 moths, soon after they 
are identified (i.e., at an apparent speed ~ 0.5c— 0.9c). 

Like in the case of PM models, in the wake of the main 
component of model OP a number of bright spots seem 
to emerge with increasingly larger apparent velocities as 
they pop up far away from the jet nozzle. However, look- 
ing at the locations from where these components seem 
to emerge, we notice that they are in clear association 
with the locus of the standing shocks of the OP mod- 
els. Such an association is even more evident when we 
look at the world-lines of the brighter features trailing 
the main component as they are localized in the 43 GHz 
radio maps. The physical origin of these trailing fea- 
tures differs from that of the trailing components seen in 
PM models. There trailing components are local incre- 
ments of the pressure and of the rest-mass density of the 
flow produced by the linear growth of KH modes in the 
beam, generated by the passage of the main hydrody- 
namic perturbation. In the beam of OP models, intrin- 
sically non-linear standing shocks are present. Nonethe- 
less, the interaction of a non-linear hydrodynamic per- 
turbation with non-linear cross shocks yields an obser- 
vational trace which resembles much that of a trailing 
component. Thereby we k eep calling such feat ures trail- 
ing components, following I Agudo et al.l (|2001f ). 

* Trailing components are pattern motions in the jet beam. 
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Fig. 17. — World-lines of a number o bright features observed in the unconvolvcd radio images resulting from the difference between 
the hydrodynamic models with and without an injected perturbation (like in Fig. II6I I. Squares, triangles and crosses correspond to the 
radio frequencies 43 GHz, 22 GHz and 15 GHz, respectively. The size of the symbols is proportional to the wavelength of the data they 
display. Shown are PM-L and OP-L (upper left and lower left panels, respectively) and PM-H and OP-H (upper right and lower right 
panels, respectively) models. The slope of the dashed line corresponds to the an apparent velocity equal to the speed of light. 



If the jet is not pressure matched, all the KH modes 
excited in the beam are blended with standing knots. In- 
deed, we realize that close to the jet nozzle, the locus of 
the first two bright spots is almost standing and, at larger 
distances, the subsequent knots show a clear increment of 
its pattern speed. The fist two trailing components are, 
actually, the traces of standing shocks which are dragged 
along with the main perturbation and oscillate around 
their equilibrium positions. The remaining trailing com- 
ponents move much faster and they can probably be due 
to the pattern motion of KH modes in the OP beam. 

Comparing the traces left by the passage of the main 
hydrodynamic perturbation in the PM and OP models 
(Fig. [T6)) . it turns out that the signatures of such per- 
turbation are much cleaner and numerous in PM than 
in OP models. The number of trailing components is 
smaller in OP than in PM models, and their world-lines 
are more oscillatory than in the latter case. For a lager 
magnetic field (models PM-H and OP-H; Fig. \W\ right 
panels) NTPs cool faster and radiate more energy, and 
thus, one can basically see only features happening close 
to the jet nozzle. 

The unconvolved data for both PM and OP models, 
and independently of the magnetic field strength, is com- 
patible with not having any time lag between the high 
and low frequency radiation emitted by the main com- 
ponent, i.e., the radiation at all three frequencies is co- 



spatial (Fig. [17]) ■ However, the convolved data display a 
number of positive and negative time lags which result 
from the difference in the size of the convolution beam 
at every frequency. In case of the PM models, there 
is a trend of the 43 GHz maximum emission to lie be- 
hind the corresponding maxima at 22 GHz and 15 GHz 
(Fig. [18] upper panels). Thereby, the low energy radia- 
tion from the main component is seen first, and later an 
observer detects radiation at higher frequencies. Never- 
theless, considering that the resolution of the convolved 
data is worse at smaller frequencies, the emission from 
the component is consistent with having no time-lags be- 
tween low and high frequency emission. This trend is in- 
dependent of the magnetic field strength, but it is more 
obvious for the model PM-H model (note the large sepa- 
ration between the different symbols beyond Zobs ~ 15i?f, 
in the Fig. llSl upper right panel). Therefore, any positive 
or negative time lag of radiation at different frequencies 
measured from convolved data has to be taken with care. 

For OP-L models, positive and negative time lags be- 
tween the high and low energy radiation are observed 
along the z-axis (Fig. [T8l lower left panel). Such time 
lags are smaller than for the PM-H model and, indeed, 
the data are compatible with no-time lags at all. For OP- 
H, in most cases, the high-frequency emission dots lie in 
front of the lower frequency ones (Fig. [18] lower right 
panel). But still, considering the difference in linear res- 
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olution for the location of the maxima, the radiation at 
different frequencies is almost co-spatial. 

Trailing components can only be tracked at 43 GHz 
close to the jet nozzle. Only after a certain distance, it 
is possible to see them at 22 GHz and even at 15 GHz 
(see the last two trailing world lines in each panel of 
Fig. [T5|) . The world lines of trailing components at 
22 GHz and, particularly, at 15 GHz, undergo substan- 
tial velocity changes. During some time intervals the 
convolved data shows receding trailing components at 
such frequencies. In the OP models, there are no clear 
trends, independent of the magnetic field strength, since 
it is very difficult to locate any local maxima at 15 GHz, 
and the 22 GHz data lie almost on top of the 43 GHz 
points. We note that there is a mismatch between the 
data points at different frequencies in the OP-H model 
at the first two recoUimation shocks (vertical threads at 
2obs ~ 4i?b and 7i?&). It is produced because there is a 
rather small relative difference in the emissivity of the 
perturbed and the quiescent jet models at 43 GHz until 
Zohs < lO-Rfe. In such conditions, the algorithm to detect 
local maxima in the space-time diagrams yields oscilla- 
tory results. A large mismatch between the world-lines 
of the peak intensity of trailing components at different 
frequencies also happens in other trailing features (e.g., 
the fourth and fifth threads in Fig. [TSl lower right panel). 
This mismatch does not exist in the corresponding un- 
convolved data (Fig. \17\ and, hence, we conclude it is an 
artifact of the finite size of the convolution beam at the 
observing frequencies. 

8. DISCUSSION AND CONCLUSIONS 

We have presented a new method (SPEV) to compute 
the evolution of NTPs coupled to relativistic plasmas 
under the assumption that these NTPs do not diffuse 
across the underlying hydrodynamic fluid. NTPs change 
their energy because of the variable hydrodynamic condi- 
tions in the flow and because of their synchrotron losses 
in an assumed background magnetic field. The inclu- 
sion of synchrotron losses and a transport algorithm for 
NTPs are major steps forward with respect to previ- 
ous approaches we have followed. The new method has 
been validated with another preexisting algorithm suited 
for the same purpose, but without including synchrotron 
losses and transport of NTPs (AM algorithm). The vali- 
dation process shows that the SPEV method reproduces 
the same qualitative phenomenology as outlined in the 
previous works of our group (G95, G97). The power of 
the new method in its whole blossom shows up when 
synchrotron cooling dominates the NTP evolution. 

8.0.0.1. Quiescent jet models: When synchrotron losses 
are considered, the resulting phenomenology can be split 
into two regimes: losses-do minated and adiabatic reg ime 
(following the convention of iMarscher fc Geailll985l ). In 
the losses-dominated regime, the knots displayed in the 
radio maps, which arc close to the jet nozzle, arc brighter 
than in models which do not include synchrotron cooling 
at the considered frequencies. Indeed, quiescent jet mod- 
els including radiative losses are more knotty than those 
models which do not include them. These features result 
from the conservation of the number density of NTPs. 
Since the same number of particles per unit of volume 
that initially extends from 7niin(^ = 0) to a certain up- 
per limit 7niax(i ~ 0) is confined into a narrower Lorentz 



factor interval, wherein more NTPs are efficiently emit- 
ting in the considered observational radio bands. In the 
adiabatic regime (reached relatively far away from the 
jet nozzle), the spectral changes, that the NTP popu- 
lation experiences as it is advected downstream the jet, 
of models with and without losses is qualitatively sim- 
ilar, since most of the high-energy NTPs (which evolve 
faster) have cooled down to energies where losses are neg- 
ligible. The beam of the jet in the adiabatic regime is 
dimmer at radio-frequencies than in models where syn- 
chrotron losses are not included. Our method lacks of 
a suitable scheme to account for diffusive shock accel- 
eration of NTPs. However, all shocks existing in the 
quiescent jet models are rather weak and, for practical 
purposes they can be considered as compressions in the 
flow, where an enhanced emission is obtained due to the 
local increase of density and of pressure. 

One of the main results of this work is that for the same 
background hydrodynamic jet model, dynamically negli- 
gible magnetic fields of different strengths yield substan- 
tially different observed morphologies. This introduces 
a new source of degeneracy (in addition to relativistic 
effects, such as, time delay, aberration, etc.) when infer- 
ring physical parameters out of observations of radio jets. 
For example, the difference in the observational proper- 
ties of models OP-L and OP-H (Sect. [O]) shows, that 
increasing the magnetic field strength by a factor of 10 
triggers a much faster cooling of the NTPs, resulting in 
a much shorter losses-dominated regime and shorter jets, 
despite magnetic field remaining dynamically unimpor- 
tant. Furthermore, jet models with such a large mag- 
netic field display a larger flux density contrast between 
shocked and unshocked jet regions. The reason being 
that after the losses-dominated regime, 7max is reduced 
so much that most of the NTP population is inefficiently 
radiating at the considered radio wavelengths and, only 
when the non-thermal electrons are compressed at cross 
shocks of the beam, they partly reenter into the efficiently 
radiating regime at the considered frequencies. 

8.0.0.2. Spectral inversion: In this paper we suggest 
that an inverted spectrum may also result if the lower 
limit of the NTP distribution 7niin is larger than the value 
of 7„ for which the synchrotron function R{x) reaches 
its maximum (Ea. 1441), in agreement w ith the theoret- 
ical predictions of iTsang fc KirkI (|2007ft . Evidences for 
flat, optically thin radio spectra in several active galactic 
nuclei have been shown by, e.g., iHuehes. AUer fc AUerl 
(|1989bf ): [Meiros3 (IT991 . and lWang. et all (|1997f) . These 
authors consider different kinds of Fermi-like acceleration 
schemes to be respon sible for the hardness of th e elec- 
tron energy spectra. IStawarz fc PetrosianI (j2008f ) show 
that stochastic interactions of radiating ultrarclativistic 
electrons with turbulence characterized by a power-law 
spectrum naturally result in a very hard (actually in- 
verted) electron energy distribution which yields a syn- 
chrotron emissivity at low frequencies with an spectral 
index ~ 1/3. Alternatively, Birk, Crusius-Watzel & 
Lesch (2001) argue that optically thin synchrotron emis- 
sion due to hard electron spectra produced in magnetic 
reconnection regions may explain the origin of flat or 
even inverted spectrum radio sources. In contrast to our 
findings, these authors, explain the spectral inversion in 
some sources as a result of a flatter electron energy dis- 
tribution. Observationally, it could be possible to dis- 
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Fig. 18. — Same as Fig. 1171 but using convolved data at every frequency. The location of local maxima trailing the main perturbation 
in the beam of PM models is difficult at 15 GHz and 22 GHz because of the large convolution beam at these two frequencies. The low 
observing resolution at such frequencies drives spurious detections of local maxima in the radio maps (trailing components), which explain 
the anomalous data at 15 GHz and 22 GHz in the range (tobsj Zobs) = (2 yr — 3.7yr, 10i?b — 20Rb) in the PM model (upper panels). The 
OP models do not display such obvious anomalies because the high intensity threads that trail the main perturbation are associated with 
preexisting recoUimation shocks in the beam, whose cmissivity, relative to the background jet, is much larger than that corresponding to 
the trailing components in the PM models. 



criminate between both possibilities by looking at the 
high-energy spectrum of the source. If there are external 
seed photons (e.g., from the AGN), which were Comp- 
ton up-scattered by the non-thermal electrons of the jet, 
the spectral index at high energies could discriminate be- 
tween the alternative explanations for the optically-thin 
inverted spectra at radio frequencies. 

Since 7inin is fixed in our model through Eq. ([55]) and 
it is not derived from first principles, one may question 
whether the value we obtain for 7min could be too large 
and, therefore, the spectral inversion we are explaining 
on the basis of taking 7min > 7m is unlikely to happen 
in nature. This would be the case if the jet was com- 
posed by an d electron-positron pla sma, in which case 
7min — 1 (^e.g. lMarscher. et al.ll2007D. F o r plas mas made 
out of electrons and protons. IWardld (|l977t ) obtained 
that for synchrotron sources with a brightness tempera- 
ture ~ 10^^ K and g = 2, 7min > 161 in order to account 
for the low degree of depolari zation in parsec- s cale e mis- 
sion regions. More recently, iBlundell. et al.l ()2006D in- 
ferr ed 7inin ^ at the hot-spots of 6C 0905-1-3955 (see 
also lTsang fc Kir3 120071 and references therein). Thus, 
the exact value of 7min is probably source dependent, and 
our minimum Lorentz factor threshold (7min — 330) can 
be well accounted by present day theory and observations 



if the jet is not a pure electron-positron plasma. 

8.0.0.3. Radio components: We have applied the SPEV 
method to calculate the spectral evolution of supcrlumi- 
nal components in relativistic, parsec-scale jets. These 
components are set up as hydrodynamic perturbations 
at the jet nozzle. For a small value of the magnetic field 
(the same as in G97), synchrotron losses are negligible 
and we recover the phenomenology shown by G97 and 
lAgudo et al] (|200ll ). 

The main component is characterized by a hardening of 
the spectrum. Pressure matched models yield a generic 
spectral profile of the component, which is rather inde- 
pendent of synchrotron losses. The hydrodynamic per- 
turbation looks in the radio maps like a burst at every 
radio-frequency and, just behind it, there is a decrease 
of the flux density. The shape of the burst is asymmetric 
in the axial jet direction, being brighter upstream than 
downstream. The shape of the burst is also frequency 
dependent because the convolution beam grows linearly 
with the observing wavelength (at lower frequencies the 
component is more symmetric in the axial jet direction). 
This triggers a decrease of the spectral index in the for- 
ward region of the main component, until it reaches a 
minimum (which precedes the intensity maxima at the 
highest observing frequency). 
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When radiative losses are important, a number of dif- 
ferences can be observed: 

1. Main component splitting in OP-H model: The 
main component splits in the radio-maps much 
more clearly than in OP-L model (Sect. [771]) . and 
the splitting takes place farther away from the noz- 
zle in the latter than in the former case. The rear 
part of the component is brighter than the forward 
one if losses are included. The spectral index pro- 
file is unaffected by the apparent splitting of the 
component. We conclude that the apparent split- 
ting of the main component is an artifact of the 
sampling of the results in the observer frame. It is 
necessary to perform a finer time sampling of the 
radio jet than the 3.5 months we have consid- 
ered in the radio maps, in which case, the main 
component exhibits an intermittent variation of its 
flux density (see on line material). If observations 
do not have the sufficient time resolution, there is 
another hint that can help to disentangle whether 
the splitting is apparent or real. In a true splitting 
of the component, each part may show a different 
spectral aging due to their different hydrodynamic 
evolutions. 

2. Radio features: Main and trailing components dis- 
play a less elongated aspect in radio maps. The 
reason is that without losses, the beam itself is 
brighter at longer distances. Thus, in the un- 
convolved data, the parts located downstream the 
jet weight more in the convolution beam than in 
the case where synchrotron losses are included. 
For the same reason, models which include syn- 
chrotron losses display a more knotty morphology. 
In the losses-dominated regime, SPEV models ex- 
hibit a brighter main superluminal component than 
SPEV-NL models. This behavior reverses in the 
adiabatic regime. Also, the ratio between the peak 
specific intensity of a trailing component to the 
specific intensity of the region of the beam imme- 
diately behind it, is larger than if losses are not 
included. The conservation of the NTPs number 
density explains this phenomenology (Sect. [5?2TT|) . 

3. Spectral properties: Behind the main component, 
the spectral index returns almost monotonically to 
its unperturbed value. In contrast, when losses are 
negligible, there is a softening of the spectrum, just 
behind the main component (where the spectral 
index reaches a maximum). 

8.0.0.4. Time lags: In this paper we explicitly show, 
that the convolved data has to be interpreted carefully. 
During most of the time the main component is observ- 
able, the radiation emitted by the component at low en- 
ergy (15 GHz) arrives to the observer before that at high 
energy (43 GHz). Indeed, for models PM-H and OP-H, a 
substantial mismatch between the world- lines of the peak 
intensity of trailing components at different frequencies 
is possible. This mismatch is an artifact due to the fi- 
nite size of the convolution beam at the observing wave- 
lengths. In contrast, the imconvolved data is consistent 
with a simultaneous emission of radiation at the three 
frequencies under consideration. This behavior matches 



our expectations, since the interval of observing wave- 
lengths is too narrow to display a substantial frequency 
dependent separation of the regions of maximum emis- 
sion. 

8.0.0.5. On the nature of trailing components: The jour- 
ney of the main component downstream the jet generates 
a number of frequency dependent bright spots which pop 
up in its wake. They differentiate themselves from the 
main component because (1) they do not emerge from 
the jet core, (2) they posses substantially s maller (some- 
times subluminal or even, receding) speeds (jAgudo et al.l 
l20nih and, as we demonstrate here, (3) they do not ex- 
hibit an obvious change in the spectral index with re- 
spect to the quiescent jet model, but (4) their observa- 
tional imprint is frequency dependent (they are clearly 
visible at the highest radio-observing frequencies, but at 
22 GHz and, particularly, at 15 GHz they are wiped out 
by the large convolution beams at this wavelengths). In 
pressure matched jet models, trailing components result 
from the linear growth of KH modes in the beam, af- 
ter the passage of t he main hydrodynamic perturbation 
(|Agudo et al.l l2001[ ). Here, we also consider overpres- 
sured jet models, where the situation is qualitatively dif- 
ferent from pressure matched ones, since the beam of 
such models develops standing shocks (non-linear struc- 
tures). Nonetheless, the interaction of a non-linear hy- 
drodynamic perturbation with non-linear cross shocks 
yields an observational trace which resembles that of a 
tr ailing component. T herefore, sticking to the definition 
of lAgudo et all (|200lh . we also call trailing components 
to the bright spots following the main component in 
ovcrpressured models, although the dynamical origin of 
such components differs. In this sense, every bright spot 
that results from the interaction between a strong hy- 
drodynamic perturbation with a relativistic beam, which 
moves slower than the main component and is not ejected 
from the jet core shall be considered as a trailing compo- 
nent. We shall add an obvious cautionary note: striving 
for the knowledge of the jet parameters, on the basis of 
a fit of the intensity variations behind a main pertur- 
bation to a number or KH modes, requires that the jet 
is pressure matched (if the jet is not pressure matched, 
all the KH modes excited in the beam are blended with 
standing knots and the predicted jet parameters might 
be inaccurate). Furthermore, it is necessary that the lin- 
ear resolution of the convolved (observational) data was 
rather good. We have tested that the unconvolved re- 
sults are roughly recovered if the FWHM of the beam at 
43 GHz is smaller than 0.25i?6. Insufficient linear reso- 
lution biases the observed features in hardly predictable 
ways, rendering inadequate the identification of features 
in the radio maps with hydrodynamic structures. 

In the future we plan to apply the SPEV method to 
perform additional parametric studies of relativistic par- 
sec scale jets. Among the parameters which can be inter- 
esting to look at, we give preference to the electron spec- 
tral index. Also, the SPEV algorithm can be coupled to 
relativistic magnetohydrodynamic codes. This will drop 
any assumption about the topology and strength of the 
magnetic field in the jet, and it will enable us to perform 
also parametric studies of polarization of the jet emission 
and of superluminal components. 
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APPENDIX 

A. IMAGING ALGORITHM 

Equations given in Sees. [3] and 3] are, in principle, sufficient to compute the synchrotron emissivity at any position 
in space and at any instant of time in the observer's frame, either using SPEV or AM methods, accounting for 
the appropriate transformations from the frame comoving with the fluid (where the emissivity [Eq. I39| , absorption 
coefficient [Eq. |42], number density of NTPs [Eq. [5], etc. are computed). The purpose of this Appendix is to explain 
the algorithm used to produce synthetic radio maps from discrete spatial and temporal elements. 

Geometry and arrival time 

While in our simulations the hydrodynamic state of the fluid is axisymmetric regardless of the jet viewing angle, the 
observed emission is, in general, not axisymmetric. We introduce the azimuthal angle (measured in the xy-plane 
from the x-axis) and define the laboratory frame (attached to the center of the AGN) 3D Cartesian coordinate system 
{x,y,z) := (i?cos0, i?sin0, Z), where the z-axis coincides with the jet axis. We denote the jet viewing angle by d, 
and choose the following observer coordinate system (rotated with respect to the 3D Cartesian system by an angle 6 
around the y-axis) 

(a;obs,2;obs, Zobs) := (x cos 6* -I- zsinfl, y, -xsin^ + zcos^) (Al) 

in which the observer is located along the Zobs axis, far from the jet. For a given elapsed simulation time T in the jet 
frame the time of observation <obs is defined as 

tohs :== T - Zobs/c (A2) 

The task of the imaging algorithm is to produce image in the (xobsi^obs) plane for a fixed arrival time tobs (note 
that the image will be symmetric with respect to the Xobs-axis if the magnetic field is completely random). From 
Eqs. lAll - IA2I it is clear that we need to have information about states of the jet at multiple instants of laboratory 
frame time in order to correctly compute the contribution at a single iobs- In a numerical hydrodynamic simulation 
we only have a finite number of discrete iterations, but each iteration has an associated time step AT. In order to 
correctly take this into account, in the following we assume that the time instant tobs has a finite duration AT as well, 
and all radiation arriving between tohs — AT/2 and tobs + AT/2 is arriving precisely at tohs- 

Particle images 

Owing to the axisymmetric nature of the problem, we only follow the Lagrangian particle motion and evolution in 
two dimensions (see Sec. [3|). However, for the purposes of imaging, a three-dimensional particle distribution needs to 
be created. We assume that each particle which is injected at the jet nozzle has a radius Ar := i?b/(2A'p), where 
i?b is the beam radius and A^p number of particles per beam radius. That means that a particle in two dimensions 
correspond to a revolution annulus in the {x^y, z) coordinate system^. In principle, by knowing the particle position 
(i?p, Zp) in the 2D grid we could compute from Eqs. lAll - IA2l all combinations of (x, y, z) and, hence, all combinations 
Xohs, 2/obs and t^hs to which the particle annulus corresponds for a fixed T. In practice, we approximate every annulus 
by a series of cubes which are distributed along a circle with radius i?p, whose center is in (0, 0, Zp). The number of 
cubes, evenly distributed in the azimuthal direction, necessary for an optimal volume coverage of the annulus depends 
on the relation between Rp and the particle radius Ar (see next subsection). By virtue of the symmetry of the jet, as 
seen by the observer, with respect to the Xobs-axis, we only need to compute the contribution from one half-annulus, 
i.e. for those cubes where y = j/obs > 0. 

We assume that both the emissivity and the absorption coefficient are homogeneous within each cube. Thus, knowing 
the particle velocity and the azimuthal angle of a given cube, we can transform its emissivity and its absorption into 
the observer frame. 

Approximation of annuli by cubes 

Given a particle with radius At and cylindrical coordinates (i?p,Zp), we approximate the corresponding half- 
revolution annulus by cubes evenly tessellating a circumference centered at (0, 0, Zp). The angular separation between 

Note that also annuli are generated from the rotation of 2D cylindrical numerical cells around the jet axis and, thereby we can apply 
the same imaging procedure when we use the cell-based algorithm AM. 
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Fig. A19. — Left: Results of the volume filling convergenec test for different numbers of particles per unit of beam radius A'^p. In the 
upper panel we show the computed length of the chord through the jet as a function of the height j/obs- We fix the value of the jet radius to 

be Rb = 0.95625. The thick black line corresponds to the analytic expectation of the chord length, i.e., 2^ ^rmb ~ fobs' ^^'^ lower panel 
the relative error with respect to the analytic expectation is displayed. Right: Results of the fiux convergence test. The figure displays the 
total fiux of images at 15 GHz for the quiescent PM (circles) and OP (squares) models as a function of the number of particles per beam 
radius Np injected at the jet nozzle. Images have been produced with 4, 8, 16, 20, 24, 32 and 40 particles per beam radius. The luminosity 



of both families of models are normalized to the total flux produced by the corresponding (PM or OP) model with TVp 



: 40. 



the cubes is defined as A(j> — min (tt, 2Ar/i?p). Thus, there are — Tr/A(f> = max (1, i?p7r/(2Ar)) cubes in a 
half-annulus of volume 14ubcs ~ N^{2Ar)'^ ~ 47ri?p(Ar)^ whereas the true volume of the half-annulus is 



V = 

' HA 



[(i?p + Ar)2 - (i?p - Ar)2]7r(2Ar)/2 = 47ri?p(Ar)2 if Rp > Ar 
(i?p + Ar)27r(2Ar)/2 = (i?p + ArfirAr if Rp < Ar 



(A3) 



In the limit of i?p > Ar, the total volume of the cubes is approximating that of the half-annulus. For i?p < Ar, 
is reset to 1 and the volume for which y > is always 4(Ar)'^, which is close to the average over all possible values 
Rp < Ar of the true volume 77r(Ar)^/6. 

Radiative transfer 

To compute an image we subdivide the (xobs,2/obs) plane into rectangular pixels, and compute the contributions to 
each pixel by checking which particle cube^° intersects which pixel at the right observation time. The ratio of the 
area of intersection to the pixel area, gives a "weight" of the contribution of a particular cube to the intensity of the 
pixel. For a given T, the value of Zobs for each particle gives the distance from the observer, so that we create a 
"linc-of-sight" (LoS) for each pixel and sort along this line all contributing particles according to Zobs (note that these 
contributions generally come from different instants of the laboratory frame time T). Since in every pixel we sum up the 
contributions spanning the observer time range [iobs — AT/2, tobs + AT/2], the intersections of every LoS with particle 
cubes are segments, not points (which would be the case if in every pixel we would only consider the instantaneous 
contributions at tobs)- After all the contributions (i.e., intersection segments) to a pixel have been accounted for, 
we solve the standard radiative transfer equation to evaluate the final pixel intensity. The above procedure can be 
performed simultaneously for a number of different values of tobs, so that a "movie" in the observer frame can be 
created. In order to transform the intensity detected in a pixel into a flux we need to multiply by the pixel area. 

Tests of the method 

In order to validate our imaging algorithm we have developed two tests which are based upon the idea that, increasing 
the number of Lagrangian particles, both the volume filling factor^^ and the total detected flux should converge. We 
first show the convergence of the volume filling method. Then we show that the images and the total flux of the 
quiescent PM-L and OP-L models converge with increasing number of particle families. 

Volume filling 

We have created a toy model consisting of a cylindrical jet with uniform velocity parallel to the jet axis, and with a 
length equal to the particle size Ar. The half-volume of such jet is Vj.i/2 = 7ri?^Ar. We inject A'p particles in the jet 

One might also use spheres instead of cubes, but we use cubes to avoid dealing with trigonometric functions and square roots when 
checking for the intersection between rectangular pixels and particles. 

We define the volume filling factor as the fraction of the jet volume occupied by our finite size Lagrangian particles. 
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evenly distributed across the jet radius (i.e., Ar = i?b/(2A^p), or _Rb = (2i + l)Ar, i = 0, . . . , Np — 1). If particles do 
not overlap, the volume filling factor is 

Ef^o' Kubc. , EtV' M2^ + l)(Ar)3 _ 1 

Since we have a finite number of particles, the jet volume is only partially patched by the volume occupied by such 
Lagrangian particles, i.e., the volume filling factor is smaller than one. Increasing the number of particles brings it 
closer to one. To test the volume filling method, we produce an "image" of the jet at an observer time tobs — with 
a 90° viewing angle, accumulating in each pixel the contributions corresponding to a laboratory frame time interval 
AT = ITbajc. However, instead of summing up the emissivity, we add up the length of the intersection of each 
particle's volume with each pixel in the (a;obs, 2/obs) plane (as described above). The idea behind the substitution of the 
emissivity by the intersection length is that, at 90° the intersection length and the intersection volume of the particles 
are proportional and, thus, measuring lengths or volumes is equivalent. 

Since we accumulate in every pixel all contributions in the range [—AT/2, AT/2], the intersection length with each 
particle equals the size of the particle perpendicular to the LoS (2Ar). Hence, the value accumulated in a pixel 
V := (a;obs,2/obs), namely Lp^, is 

^P- = ET^2Ar, (A5) 

where Ai and are the area of intersection of a particle with a pixel and the pixel area, respectively. The sum in 
Eq. (jA5p extends over all particles that are intersected by the line of sight that departs from P. In the limit Ar 
(equivalently, TVp oo) Ai 4(Ar)^. On the other hand, the number of particles intersected by the LoS departing 
from V and having a cross sectional area Ap^ is Np^ = Apx/(2Ar)^. Therefore, we have 

i-„^P^ = i-o«^^^)'£ - ^^i^l-vlJ . (A6) 

Equation [X6] simply expresses that, in the limit A^p — s- oo, the length measured in the pixel V should tend to the length 
of the chord determined by the intersection of the jet body with the line of sight from V. Figure [Al9b shows that for 
Np > 16 the results converge very rapidly to the analytic expectation (thick black line). 

Total flux 

To test the convergence of the imaging algorithm we have produced images of quiescent PM-L and OP-L models 
with varying Np. The total number of particles in the grid grows as Np, it is thus important to minimize the number 
of particle families for numerical purposes. On Fig. IA19b we show the total image flux at 15 GHz for PM-L and 
OP-L models as a function of A'p. The values are normalized to the flux of the model with the largest number of 
injected particles per beam radius (A^p = 40), which we consider the reference value. This test is important because 
the total flux represents a global value of every model, since it is computed by summing up the individual fluxes 
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Fig. A21.— Same as Fig. lX20l but for the OP-L model. 



arriving to each pixel in the detector, and multiplying by the corresponding pixel area. Remarkably, for A'p > 16 
the flux does not deviate more than 5% form the reference value. Thus, any model with TVp > 16 has sufficiently 
converged to an appropriate total flux. This has motivated our choice to work with iVp — 32 in the current paper, 
since it yields an optimal trade-ofF between numerical accuracy and computational cost. Figures IA20I and IA21I show 
images corresponding to the convergence tests for models PM-L and OP-L, respectively. 
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